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Abstract 

This thesis presents a novel design methodology for the synthesis of automatic 
controllers, together with a computational environment — the Control Engineer's 
Workbench — integrating a suite of programs that automatically analyze and de- 
sign controllers for high-performance, global control of nonlinear systems. This 
work demonstrates that difficult control synthesis tasks can be automated, using 
programs that actively exploit and efficiently represent knowledge of nonlinear dy- 
namics and phase space and effectively use the representation to guide and perform 
the control design. The Control Engineer's Workbench combines powerful numer- 
ical and symbolic computations with spatial reasoning techniques. The two major 
programs in the Workbench — Phase Space Navigator and MAPS — work together 
to model and reason about the phase-space geometry and topology of a given sys- 
tem, to plan global control reference trajectories, and to navigate the system along 
the planned trajectories. They use a novel technique of "flow pipes" to group in- 
finite numbers of distinct behaviors into a manageable discrete set that becomes 
the basis for establishing the reference trajectories. 

As a demonstration of this approach, I exhibit the automatic design of a non- 
linear controller for a magnetic levitation system. The control system synthesized 
by the Workbench can stabilize a maglev vehicle with large initial displacements 
from an equilibrium and outperform the classical linear feedback design for the 
same system by a factor of 20. 
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Chapter 1 



Introduction 



The purpose of computing is insight, not numbers. 

— R. W. Hamming 

Synthesis and analysis of control systems for complex physical systems are 
difficult tasks. Complex systems operate in large nonlinear regimes; they often 
comprise many components and contain a large number of state variables. Linear 
systems have been well studied, and there are systematic design procedures for lin- 
ear control systems. On the other hand, nonlinear systems are far less understood. 
Nonlinear control synthesis is limited by the lack of available control schemes and 
analysis techniques; traditional design and analysis are seldom applicable to these 
systems. 

Powerful computers can help engineers and scientists unveil the intricacies of 
nonlinear phenomena and can revolutionize the design of complex, high-performance 
control systems. The computational exploration of nonlinear systems has drasti- 
cally changed the way we think about the world. This opens the possibility of 
simulating and designing physical systems, without resorting to linear approxima- 
tions. 

Yet despite advances in computer technology, nonlinear control systems are 
still difficult to design and analyze. The difficulties arise from the lack of design 
methodologies that actively exploit the special nature of nonlinearities, and from 
the lack of high-level computational tools that can harness the available compu- 
tational power. Few nonlinear systems admit closed-form solutions. The ability 
to understand and control these systems requires extensive numerical simulations. 
The burden of the task is on human engineers to carefully design experiments, to 
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distill qualitative information from numerical simulation results, and to use the 
information to design appropriate control actions. This manual, numerical ap- 
proach is burdensome and error-prone, and severely limits the design space one 
can explore. 

The modern geometric theory of dynamical systems, pioneered by Poincare 
at the beginning of the century, provides a qualitative way to describe the rich 
dynamical behaviors of nonlinear systems. The geometric representation forms 
a viable substrate for computationally analyzing qualitative aspects of dynamics 
and exploring novel control-synthesis strategies. Programs that employ powerful 
symbolic and numerical techniques have already shown promise in the geometric 
analysis of nonlinear dynamical systems [2]. 

However, the phase-space geometry and topology of a dynamical system are 
hard to describe in a way that allows for efficient computational manipulations for 
the purpose of control design. Nonlinear systems can have extremely convoluted 
phase-space geometries; the complexity becomes much worse as the dimensional- 
ity increases. Humans can comfortably picture and manipulate two and three- 
dimensional objects with the aid of graphic, geometric modeling techniques. For 
higher-dimensional systems, however, few visualization and manipulation tools ex- 
ist. Automatic modeling, analysis, and design tools are necessary to identify, 
extract, and reason about the spatial properties of phase space. 

1.1 Contributions 

This thesis demonstrates that difficult control synthesis tasks can be automated, 
using a computational workbench consisting of programs that actively exploit 
knowledge of nonlinear dynamics and phase space. These programs combine pow- 
erful numerical and symbolic computations with spatial reasoning techniques. The 
thesis contributes to the state of the art in artificial intelligence and control systems 
engineering in several ways: 

• The thesis has developed a qualitative representation for complex behav- 
iors of dynamical systems in phase space and a design language for com- 
putationally expressing and manipulating these behaviors. The qualitative 
representation captures the gross aspects of dynamics in a relational graph 
of phase-space structure and a set of discrete objects called flow pipes — the 



equivalence classes of behaviors. The design language describes a control de- 
sign task in terms of well-defined geometric, combinatorial operations on the 
flow pipes. This language helps formalize aspects of implicit expert reasoning 
of control engineers in solving control design problems. The representation 
and the language are developed independently of the orders of systems, i.e., 
the dimensionality of phase spaces. 

The thesis has constructed a computational environment, the Control En- 
gineer's Workbench, integrating a suite of programs that automatically 
analyze and design high-performance, global controllers for a large class of 
nonlinear systems. These programs combine powerful techniques from nu- 
merical and symbolic computations with novel representation and reason- 
ing mechanisms of artificial intelligence. The two major components in the 
Workbench — the Phase Space Navigator and MAPS 1 — work together to vi- 
sualize and model the phase-space geometry and topology of a given system. 
They reason about and manipulate the phase-space geometry and topology 
and search for optimal control paths connecting initial state and the desired 
state for the system. The Workbench represents the result of design and 
analysis in a symbolic form manipulable by other programs, and produces 
a high-level summary meaningful to professional engineers. It also presents 
the result in a graphical form. 

The thesis presents a novel phase-space design methodology for the syn- 
thesis of automatic controllers. The design methodology computationally 
exploits dynamical systems' nonlinearities in terms of phase-space geome- 
tries and topology. It designs a prespecified control law — control reference 
trajectories — for a system by synthesizing the desired shape for phase-space 
geometry dictating trajectory flows. It uses the novel technique of flow pipes 
to group infinite numbers of distinct behaviors into a manageable discrete set 
that becomes the basis for establishing reference trajectories, and navigates 
the system along the planned reference trajectories. The phase-space design 
approach requires powerful computational tools that are able to identify, ex- 
tract, represent, and manipulate qualitative features of phase space, and is 
embodied in programs comprised in the Workbench. 



1 MAPS stands for Modeler and Analyzer for Phase Spaces. 



• The thesis has demonstrated the Workbench and the phase-space design 
methodology in an application of great practical interest: the Workbench 
helped design a high-quality controller for a magnetic levitation system — the 
German Transrapid system. The control system synthesized by the Work- 
bench stabilizes a maglev vehicle with much larger initial displacements than 
those allowed in a previous linear design for the same system using classical 
linear feedback technique. The simulation shows that our design outperforms 
the linear design by a factor of 20. Professional control engineers consider 
this result important enough to present at the 31st IEEE Conference on 
Decision and Control [56]. 

1.2 Goal: Mechanizing Control Design 

A real-world control system is a complex closed-loop system with extremely rich 
dynamics. Computation and reasoning are pervasive in the design and operation 
of the controller. Sensors collect a large amount of quantitative information. State 
and parameter estimators infer hidden information about the system from the 
sensed data. The system is modeled with a representation appropriate for further 
analysis and design based on available information. The model is then analyzed 
to extract behaviors that are judged significant for the control objective. To meet 
the control objective, a control law is synthesized to change the natural dynamics 
of the system. 

The domain of automatic control brings together issues of sensing, estimation, 
control synthesis, and control execution. The study of their common themes — 
computation and reasoning — serves as a framework for coherently addressing 
these issues and makes it possible to employ advanced computational techniques 
to drastically improve modern control design. This thesis focuses on the control 
synthesis. 

A control engineer goes through the following design steps to synthesize a 
controller for a given physical system: 

1. Modeling: translate the physics and the constraints of the physical system 
into a quantitative mathematical model. 

2. Analysis: analyze the model of the system. 

3. Design: design a controller for the system. 



4. Verification: verify the control design. Iterate Steps 2-4 if necessary. 

5. Implementation: implement the control design in software or hardware. 

In the above control design procedure, Step 1 models the system, usually with a 
set of differential equations. Step 2 examines the model and analyzes the behaviors 
of the system. Based on the analysis, Step 3 arrives at a control design according 
to the prespecified control constraints. Step 4 ensures that the design meets the 
control specification. The last step implements the control design to control the real 
physical system. The Steps 2, 3, and 4 are often iterated before a reasonable control 
law is synthesized. Except for very few cases in which analytic-form solutions are 
available, computer simulations are the main tools for analyzing the systems and 
for designing and verifying the controllers. 

The task of control design and analysis requires powerful knowledge rep- 
resentation and reasoning mechanisms and computer simulation tech- 
niques. The mechanization of this task and the exploration of novel 
control design methods are the subjects of this thesis research. 

1.3 Limitations of Traditional Methods 

We are interested in developing a computational design method suitable for dealing 
with the complexities of real-world nonlinear control systems and in mechanizing 
the design process by computer programs. Traditional methods relying on purely 
analytic or numeric techniques or on linear design theory are inadequate for car- 
rying out and automating the design process. 

Systematic techniques for control design of linear systems have been well de- 
veloped. Controllers for linear systems can be synthesized via methods such as 
pole-placement or Bode diagrams [35], with satisfactory performance. However, 
complex control systems operate in larger regions where system nonlinearities are 
too significant to ignore. These nonlinear control systems do not lend themselves 
easily to linear (small signal) analysis and design. There have been attempts to 
extend linear control techniques to nonlinear control systems: a piecewise-linear 
control design, also called gain scheduling, approximates a nonlinear function with 
several linear pieces, each of which admits a linear control design. For example, 
in aircraft autopilot design, an airplane system model is linearized around several 



hundred operating points selected within the plane's operating region. Each oper- 
ating point is subject to a linear controller design. However, this type of control 
can be very expensive and complex in both design and implementation. For a 
highly nonlinear system operating in a large region, an enormous number of lin- 
earizations and controller designs have to be carried out. The complexity and the 
cost of conventional implementations prohibit many practical applications. 

The computational method developed in this thesis designs a nonlinear control 
system in phase space. Unlike gain scheduling, our method explicitly models and 
actively exploits the nonlinearities in terms of trajectory flows and can carry out 
control design automatically. Although phase space is an extremely useful medium 
for exploring novel nonlinear design methods, the phase-space design approach 
requires sophisticated computational techniques and representational mechanisms 
to make the approach computationally feasible. 

Existing control simulation softwares are inadequate for automatically design- 
ing highly complex nonlinear systems. Commercially available programs like MAT- 
LAB and SlMULAB [30] rely on numerical simulations. These programs are, at their 
very best, semi-automatic and serve as interactive design aids to human engineers. 
Although they are equipped with elaborate graphic interfaces, these programs 
provide only fragmented, limited capabilities such as integrations and equilibrium 
determination for performing the simulation task; human users need to prepare the 
simulation and to interpret the result. The specialized control toolboxes embedded 
in these programs are "shallow" expert systems; they lack deep domain knowledge 
and do not have mechanisms for computationally representing and manipulating 
a control design. 

1.4 Characteristics of the Approach of This Thesis 

The computational approach presented in this thesis has the following character- 
istics: it can analyze and design nonlinear, global control laws; it is autonomous; 
and its development is independent of the order of a system. 

Our approach exploits the richness of nonlinear dynamics and phase-space rep- 
resentations. For example, the programs use the flow pipes to explore a variety 
of control paths otherwise difficult to obtain with just a piecewise linearization. 
Our approach designs control trajectories from a global point of view and thereby 
maximizes the effect of control actions to obtain good performance. 



The programs embody significant knowledge of dynamical systems theory and 
control theory. Because of their autonomous nature, the programs can perform 
certain tasks that human counterparts have difficulty with, for instance, the task 
of visualizing and manipulating complicated high-dimensional geometric objects. 
The machinery for geometric modeling is developed independently of the spatial 
dimensionality. Hence the programs are not restricted to low-order systems, as 
long as the computational complexity allows. 

These programs are mostly useful in synthesizing high-performance control 
systems that do not lend themselves to traditional design and analysis techniques. 
On the other hand, the qualitative analysis could also assist engineers in exploring 
novel control design methods. Aside from the utility of these programs, this thesis 
formulates the task of control design in a computational way; this formulation 
provides a concrete substrate and a rich domain for exploring and developing new 
reasoning and representation techniques otherwise hard to grasp in traditional 
computer science and artificial intelligence research. 

1.5 What the Control Engineer's Workbench Does 

Given a model of a physical system and a control objective, the Control Engineer's 
Workbench analyzes the system and designs a control law achieving the control 
objective. A user typically interacts with the Workbench in the following way. 

The user first tells the Workbench about the system: he inputs a system model 
in terms of a differential equation, parameter values, and bounds on state variables 
for analysis in the form of a phase-space region. The user also tells the Workbench 
about the requirements on the control design: he specifies the desired state for the 
system to settle in, the initial states of the system, the allowable control parameter 
values, and the constraints on the control responses. 

The user then asks the Workbench to analyze the system within the parameter 
ranges of the model. The Workbench visualizes the totality of the behaviors of 
the system over the parameter ranges; it represents the qualitative aspects of the 
system in a data structure and reports to the user a high-level, symbolic summary 
of the system behaviors and, if necessary, a graphic visualization of the phase-space 
qualitative features. 

Next, the user instructs the Workbench to synthesize a control law for the sys- 
tem, subject to the specified design requirements. The Workbench searches for 



the global control paths that connect the initial states of the system and the de- 
sired state, using the qualitative description about the system. More specifically, 
the search is conducted in a collection of discrete entities representing trajectory 
flows in phase space. After the global control paths are established, the Work- 
bench determines the controllable region of the system and the switching surfaces 
where control parameters should change values. A synthesized control reference 
trajectory consists of a sequence of trajectory segments, each of which is under 
a constant control. For a point-to-point control design, the Workbench further 
constructs a smoothed trajectory connecting an initial state and the desired state. 

The Workbench has designed a high-performance nonlinear control law for an 
electro-mechanical system of practical interest — the stabilization of the German 
Transrapid magnetic levitation transportation system. The details of this design 
will be described in Chapter 5. Here, we briefly illustrate what the Workbench 
does in this application. 

The train and guideway system of the Transrapid is described by a third-order 
nonlinear differential equation: 

dx z(Vj—Rx) _|_ an/ 

dt LqZq z 

dy — n _ LqZqx 2 
dt & 2mz 2 

dt y 

Ensuring smooth ride-quality for passengers is one of the top priorities in the 
design of the Transrapid system. Unfortunately, the system is unstable at the 
operating equilibrium for parameter V{ at 140 volts. The objective of the control 
design is to synthesize a control law that stabilizes the system at the operating 
equilibrium. 

However, the stabilization control of the train, typically traveling at about 350 
miles per hour along the guideway subject to various disturbances, is a very difficult 
task, for the system operates in a large operating region where nonlinearities come 
from the geometries of the train and the guideway and from the magnetic forces. 
The system requires a high-quality global stabilization design. 

The Workbench is given the system equation in a symbolic form, the ranges 
of parameter values, a bounding box in phase space, and the control objective. 
It analyzes the behaviors of the system for the given parameter values and syn- 
thesizes a controllable region about the previously unstable operating equilibrium; 
the region is maximum with respect to the allowable control authority. 
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(a) 




(b) 



<equilibrium-points : 
equilibrium 1. (saddle at (140. 0. 200. ))> 

<f low-pipes : 
flow-pipe 1. (from *infinity* to *infinity*) 

boundary: (2D-manif old) 
flow-pipe 2. (from *infinity* to *infinity*) 

boundary: (2D-manif old)> 

<controllable-region : 
interior: 831 tetrahedra 
boundary: 526 triangles 
volume: 0.848 

range-of -control: (140 300) 
number-of -switching-surf aces: 1 
maximum-displacement-z: 0.00455> 



(c) 



Figure 1.1: The high-level description of the control design for a maglev system, 
automatically generated by the Workbench: (a) a projection of the controllable 
region showing several synthesized reference trajectories and the operating equilib- 
rium; (b) a graphical rendering of the controllable region; (c) a symbolic summary 

characterizing the system and the control design. 

9 



The Workbench reports the results of the analysis and the control design in a 
graphic and symbolic description in Figure 1.1. The description contains graphic 
displays of the controllable region showing the region boundary, the equilibrium, 
and synthesized control reference trajectories from sampled initial states; it also 
characterizes the behaviors of the system and the controllable region in a symbolic 
form. Although the pictures shown here are in black-and-white, the actual screen- 
display for the region produced by the Workbench is in color. The result of this 
application has been presented to professional control engineers and been favorably 
compared with a previous manual, linear design for the system [56]. 

1.6 Organization of the Thesis 

The rest of the thesis is organized as follows. 

Chapter 2 first gives an overview of the Control Engineer's Workbench. The 
bulk of the thesis is then divided into three parts: Chapter 3 develops an analysis 
method, MAPS, for constructing qualitative representations of dynamics in terms 
of asymptotic behaviors and equivalence classes, i.e., the flow pipes; Chapter 4 
presents Phase Space Navigator for the synthesis of control laws in phase space. 
The Phase Space Navigator searches for optimal control trajectories in a collec- 
tion of flow pipes forming a flow-pipe graph. Chapter 5 applies the machineries 
developed in Chapters 3 and 4 to an engineering problem of practical interest: the 
design of a global, nonlinear controller for the stabilization of a magnetic levita- 
tion transportation system. Finally, Chapter 6 describes other related work in the 
computational analysis and design for control systems, and Chapter 7 concludes 
the thesis with a summary and future work in the area. 

Although the domain this thesis studies is control design, readers with knowl- 
edge of linear algebra, elementary differential equations, and an introductory to 
linear systems or control will be able to go through the thesis with little difficulty. 
I have included necessary background materials on dynamical systems theory and 
control theory to make the thesis self-contained. 
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Chapter 2 

The Control Engineer's 
Workbench 



This chapter gives an overview of the Control Engineer's Workbench — an imple- 
mented computational environment for analyzing and designing dynamical control 
systems in phase space. The Workbench, consisting of a suite of programs, allows 
engineers to visualize the phase-space structures of the systems, to manipulate the 
dynamics, and to synthesize controllers. 

2.1 Reasoning About Control Design 

How does a control expert reason about a control design task? A professional 
control engineer uses a body of specialized knowledge to carry out the design. 
The engineer develops insight through an analysis of the physical system, uses the 
insight to explore the design space constrained by control requirements, and makes 
engineering judgment about design choices and trade-offs. A particularly intuitive 
design method is the phase-plane method for analyzing and designing a 2nd-order 
nonlinear control system in a phase plane [35]. 

Let us consider an example of the deformation of an elastic column. Putting 
a weight on top of the column, the column oscillates around its principal axis. 
Putting a heavier weight, the column oscillates and settles to a buckled state 
after a short while. How would an engineer analyze the behavior of the system 
and design a control law to stabilize the buckling motion with the phase-plane 
method? 

The engineer first describes the system with a mathematical model. Compiling 
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the dynamics into a suitable model is a nontrivial task. Discussion of model build- 
ing is outside the scope of this thesis. We assume that the engineer is given a simple 
model for the column in the form of an equation of motion, a 2nd-order nonlinear 
ordinary differential equation. Starting with the model and a set of typical pa- 
rameter values, the engineer determines all possible responses of the system within 
certain ranges of amplitudes and input signals. The responses with different initial 
conditions from the simulation are then plotted in a two-dimensional plane — the 
phase plane (Figure 2.1). A state, a pair of position and velocity, is plotted as a 
point in the plane. A trajectory consists of a sequence of consecutive states and is 
drawn as a one-dimensional curve. With this graphical plot, the engineer identifies 
the natures and locations of critical points. He then inspects the gross features of 
trajectories: the trajectories approach two critical points in spirals, respectively. 
The engineer groups the trajectories into two classes, each of which corresponds 
to those approaching the same critical point. The phase plane gives a qualitative 
picture of the system responses. The shaded region in the picture is the collec- 
tion of trajectories going to the left-hand critical point, one of the buckled states. 
The critical point at the middle of the picture is the unbent state of the column. 
In order for the engineer to interpret the responses, he needs to know things like 
system states, trajectories, critical points, and qualitative regions. 

The buckling motion that leads the column to a buckled state traces out a 
trajectory in the phase plane, for instance, a trajectory that starts at a location in 
the shaded region and evolves to the left-hand critical point in Figure 2.1. Suppose 
the engineer wants the column to return to its original unbent state after the oscil- 
lation dies down, instead of letting it get buckled. He introduces a control input, 
a force exerting on the column, and observes how the system responses vary with 
the input. With the pictorial description of the responses, the engineer determines 
a solution curve in phase space as a control law that moves the trajectory to the 
middle critical point. 

In short, the engineer goes through the following steps with the phase-plane 
method: 

• simulate the system extensively with different initial conditions 

• plot the behaviors in trajectories on a piece of paper 

• interpret the result with visual inspection 
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Figure 2.1: The phase plane of the buckling column (from Abraham Sz Shaw 1988). 

• design a control law to obtain desired behavior. 

The phase-plane method is able to analyze the transient and asymptotic behav- 
iors of planar systems, whether linear or nonlinear. Kalman had used this method 
to design switching control laws based on the phase-plane plot [26]. His method 
decomposes a phase plane into discrete linear regions, designs linear control law 
for each region, and matches transients on the boundaries of the regions. 

The phase-plane method is a useful tool for analyzing qualitative responses of 
a control system. It is, however, manual, prohibitively expensive, and confined 
to planes. Although the diagrammatical sketch of a phase plane illustrates just 
the qualitative aspects of the systemj to obtain a picture like this requires lots of 
human effort in preparing numerical simulations, collecting the numbers, sketching 
the results, analyzing the trajectory plot, and interpreting it in a qualitative picture 
like Figure 2.1. Worse, this method becomes useless in cases when the order of a 
system is greater than two and the n^nlinearity results in convoluted phase-space 
geometry. The mechanization of the task with autonomous computer programs 
would alleviate many of these restrictions. 
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name : 


buckling_column 




equation_of -motion: 


J dx-i/dt = x-2 

1 dz'ijdt = —p\Xi - 


- p 2 x\ - p 3 x 2 + u 


state_variable : 


X\ 




state_variable : 


x 2 




parameter: 


Pi = 


-2.0 




parameter: 


P2 = 


1.0 




parameter: 


P3 = 


0.2 




parameter: 


U = 


0.0 




bounding_box : 


xi G [-3.0,3.0], 


x 2 e [-4.0,4.0] 



Figure 2.2: Engineer Input: the model for a buckling steel column. 

2.2 A Real Scenario with the Workbench 



The Control Engineer's Workbench automates a significant portion of the control 
engineer's design task. In the following session, the Workbench autonomously 
analyzes the buckling motion of a steed column under compression and synthesizes 
a control law to stabilize the motion. The interaction between a human user and 
the Workbench is annotated in explanatory italics for ease of reading. 



• The engineer types in the model of the elastic column buckling 
under axial compressive force and asks the Workbench to analyze 
the system for the given parameter values (see Figure 2.2). 

• The Workbench summarizes the qualitative behaviors of the sys- 
tem in terms of equilibria, stability regions, and flow pipes (see- 
Figure 2.3). 

The Workbejich also displays the phase portrait of the system 
showing the first flow pipe (see Figure 2-4). 

• The engineer instructs the Workbench to design a control law that 
stabilizes the buckling motion, in particular, a law that brings the 
column back to the unbuckled state corresponding to the unstable 
equilibrium at the center $f Figure 2.4 . He specifies the control 
design requirements (Figur^p 2.5). 
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< equilibrium-points : 




equilibrium 1. (attractor at (1.41 p.)) 


equilibrium 2. (saddle at (0. 0.)) 




equilibrium 3. (attractor at (-1.41 


o.))> 



<trajectories : 

<boundary-traj ectories : 
trajectory 1. (from *infinity* to (o. 0.)) 
trajectory 2. (from *infinity* to (0. 0.))> 

<connecting-traj ectories : 
trajectory 3. (from (0. 0.) to (-1,41 0.)) 
trajectory 4. (from (0. 0.) to (1.41 0.))» 



<stability-regions : 
stability-region 1. 

attractor at *infinity* 

stability-boundary: () 

connecting-trajectories: () 
stability-region 2. 

attractor at (1.41 0.) 

stability-boundary: (trajectory 2. trajectory 1.) 

connecting-trajectories: (trajectory 4) 
stability-region 3. 

attractor at (-1.41 0.) 

stability-boundary: (trajectory 2. trajectory 1.) 

connecting-trajectories: (trajectory 3)> 

<f low-pipes : 
flow-pipe 1. (from *infinity* to (-1.41 0.)) 

boundary: (trajectory 2. trajectory 1. trajectory 3.) 
flow-pipe 2. (from *infinity* to (1,41 0.)) 

boundary: (trajectory 2. trajectory 1. trajectory 4.)> 



Figure 2.3: Workbench Output; a symbolic summary of the analysis. 
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Figure 2.4: Workbench Output: the phase portrait of the buckling column showing 
equilibrium points, boundaries of stability regions, and a flow pipe leading to the 
left-hand attractor. 



control_type: 
goal_state: 
initial _state: 
range_of -control : 



point_to_point 
(0.0, 0.0) 
(-1.0, -3.0) 
u e [-0.2,0.2] 



Figure 2.5: Engineer Input: the control objective. 
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The Workbench reports a synthesized control reference trajectory 
(see Figure 2.6 for a phase-space illustration and a script of the 
control law). The global portion of the reference trajectory shown 
in the figure consists of four segments, each of which starts at a 
switching state marked as a small circle. The control parameter 
value is held constant for each segment and is denoted by U\, U2, 
U3, or Ui. The control law for the reference trajectory is rep- 
resented as a sequence of tuples: (time, switching-state, control). 
The control law invokes a local controller when the trajectory en- 
ters the vicinity of the goal. The Workbench determines that the 
control response time for the global portion of the reference trajec- 
tory is less than 7 seconds^ 



2.3 Overview of the Workbench 

In the previous session with the Workbench, the Workbench has demonstrated the 
following capabilities: 

• Deciding what behaviors are significant. The Workbench looks for qualitative 



features like equilibrium points, 



stability regions, and trajectory flows. 



• Describing the behaviors qualitatively in computational terms. The Work- 
bench models the stability regions and trajectory flows geometrically. 

• Reasoning about a phase space with flow pipes. 

• Performing the control design autonomously. 

These capabilities are supported in the Workbench by (1) qualitative represen- 
tation consisting of dimension-independent constructs, (2) hierarchical extraction 
of the behaviors, (3) modeling and manipulation mechanisms for flow pipes, and 
(4) algorithms implementing the geometric, combinatorial, and numerical compu- 
tations. 



2.3.1 Design requirements 

We motivate the design criteria for the Workbench with the intended task domain 
of control synthesis and with computational considerations. 
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; ; The Synthesized Control Law: specifying the time instance, switching 
;; state, and corresponding control Value for each switching: 



((time 0.) (switching-state #(-1 -3)) (control .2)) 

((time .284) (switching-state #(-1.82 -2.71)) (control 0.)) 

((time 1.06) (switching-state #(-1.80 2.49)) (control -.2)) 

((time 2.71) (switching-state #(1.35 1.82)) (control 0.)) 

((time 6.76) (switching-state #(-.00^3 -.0692)) (control *local-control*) ) 



Figure 2.6: Workbench Output: the control reference trajectory and the control 
law for stabilizing the column. In the t>lot, the reference trajectory is drawn in solid 
lines, and the flow-pipe boundaries of the uncontrolled column in dashed lines. 



18 



• The design of the Workbench CcJls for a concise representation that captures 
essential features of a control system. The representation should be sensible 
to human engineers and manipulable by other programs. The Workbench 
needs to present the result of its analysis and design to human designers 

to interact with the design in a direct way. 
This communication requires a high-level, intuitive presentation of the result. 
Other programs in the Workbench need to efficiently access and manipulate 
the representation. Instead of encapsulating everything about the system, 
the representation should contain only information that is useful for the 
control synthesis task. We have fchosen a qualitative representation describing 
a system in terms of its phase-sjpace geometric features for this purpose. 

• The Workbench needs modeling algorithms to efficiently construct the rep- 
resentation. The algorithms should identify and extract implicit dynami- 
cal properties from numerical explorations and summarize the result in the 
qualitative form. In the implementation, the Workbench internally uses a hi- 
erarchy of intermediate representations. The qualitative information about 
the system is extracted in a step-by-step fashion, from local descriptions to 
global ones. 

• The Workbench also needs a reasoning mechanism to efficiently manipu- 
late the representation for synthesizing a control law. The program searches 
through the representation to f|nd feasible control trajectories. The Work- 
bench uses a graph mechanism that manipulates a discrete collection of ob- 
jects called the flow pipes. 



2.3.2 Anatomy of the Workbench 

The Workbench analyzes and designs control systems in phase space. In the sce- 
nario in Section 2.2, the Workbench locates equilibria and models stability regions; 
it groups trajectories, i.e., behaviors, }nto the flow pipes; it constructs a dimension- 
independent representation for the phase space; it designs control reference tra- 
jectories with this representation. In addition, the Workbench communicates with 
the user in high-level terms through an interface for inputting the model and for 
describing the design and analysis and a graphic module for displaying the result. 
Like the phase-plane analysis, the Workbench is able to reason about and ma- 
nipulate dynamics in terms of phase-space geometries. It partitions a phase space 
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into discrete regions. For the previous elastic column example, the program au- 



tonomously explores the behaviors of 
ture of the phase portrait in Figure 2 
kind of information as that of Figure 



the system and generates the high-level pic- 
.4. The picture contains essentially the same 
2.1. 
But unlike the phase-plane method, the Workbench provides computational 
means for manipulating the dynamics. It decomposes the phase space into sub- 
regions that can even be globally nonlinear. The shaded region in Figure 2.4 is 

The program internally represents the criti- 
cal points and geometries of regions in a data structure that allows other programs 
to manipulate, visualize, and communicate with human users. Because the geome- 
try and topology of a system's phase sjpace are modeled with a simplicial structure, 
the representation and reasoning mechanisms for this structure are independent of 
the dimensionality of phase space. 

The Workbench serves as an intelligent assistant to control engineers. The 
components of the Workbench are shown in Figure 2.7: 

• the MAPS program for simulation and interpretation 

• the Phase Space Navigator for control synthesis 

• a graphic program for visualizing the design 

• a user interface for communication with the system. 

The component for model building is inot in the Workbench yet. The programs in 
the Workbench implement various algorithms: symbolic differentiation, numerical 
algorithms on differential equations, numerical integrations, modeling of geomet- 
ric structures, clustering of equivalence classes, graph algorithms, etc. The infer- 
ence mechanism of the Workbench uses these programs to construct a qualitative 
phase-space structure for representing a system, to check for the consistency of the 
structure, and to reason about and manipulate the representation through a graph 
of flow pipes. 

The flow of computation within \he Workbench is illustrated in Figure 2.8. 
Given a model of a system, a bounded phase-space region of interest, allowable 
parameter values, and control objectives and constraints, the Workbench performs 
stability and trajectory flow analysis for the system in phase space and interprets 
the result in a phase-space graph. Th(£ Workbench then explores the control space 
to synthesize a desired control law subject to the design constraints. It reports 
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The Control Engineer's Workbench 
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Figure 2.7: The structure of the Workbench. 
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a point-to-point control, the output 
consists of a sequence of tuples of t 
control value; if an initial operating 



the control law specifying reference trajectories and performance properties. The 
control design for steering towards an equilibrium is performed in this way: for 

is a reference trajectory, whose control law 
:ime, switching state and the corresponding 
region is given, the output is a controllable 
region geometrically represented as a! polyhedral structure. The steering towards 
a limit cycle has not been implemented yet, although its implementation would 
be very much the same in Poincare [sections as the case for steering towards an 
equilibrium. 

The programs in the Workbench ^re written in Scheme, a dialect of LISP. The 
Workbench is implemented on HP 70Q series workstations. For moderately difficult 
problems, for example, the buckling steel column example, it takes the Workbench 
on the order of minutes to arrive at a control law. 
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Figure 2.8: The flow of Computation in the Workbench. 
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Chapter 3 



Automatic 
Modeling and 



Phase 



-Space 
Analysis — MAPS 



This chapter develops and demonstrates a computational method for automati- 



cally analyzing qualitative behaviors 



of dynamical systems in phase space. We 



will discuss the extraction and representation of the qualitative behaviors of a dy- 
namical system in terms of a qualitative phase-space structure describing spatial 
arrangement of geometric entities in pjhase space, a geometric construct called flow 
pipe for modeling trajectory flows, aiji analysis algorithm and its implementation 
in the MAPS program, and examples illustrating the mechanism of the analysis. 

3.1 Introduction 

Analysis of dynamical systems via p|ase-space geometric structures plays an in- 
creasingly important role in experimenting, interpreting, and controlling complex 
systems [1, 2]. Nonlinear systems usjually fall outside the domain of traditional 
analysis methods, such as Fourier analysis for linear systems. However, most of 
the important qualitative behaviors of a nonlinear system can be made explicit in 
phase space with a geometric analysis [39]. 

MAPS is a program for understanding and representing qualitative structures 
of phase spaces. MAPS combines numerical, symbolic, and geometric compu- 
tations with techniques of spatial reckoning. It embodies theoretical knowledge 
about nonlinear dynamical systems and formulates in computational terms a sig- 
nificant amount of informal working knowledge of professional control engineers in 
analyzing complex control systems, in particular, in the form of the phase-plane 
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method [35]. We will illustrate our techniques of extracting and representing the 
qualitative phase-space features with two- and three-dimensional systems. The 
techniques presented in this chapter also apply to higher-dimensional dynamical 
systems. 

Complex systems are often nonlinear and high dimensional. Our theoretical 
knowledge about nonlinear dynamical systems is far from complete. Therefore, 
many engineering applications rely on extensive numerical experiments. A numer- 
ical simulation typically generates an immense amount of quantitative information 
about a complex system. To interpret the numerical result and to use the in- 
formation for engineering designs, it is essential to develop qualitative methods 
that automatically analyze the system, extract the qualitative features, and rep- 
resent them in a high-level description sensible to human beings and manipulable 
by other programs. The representation for the qualitative behaviors needs to be 
parsimonious and yet capture the essential properties of the dynamical system un- 
der study. The task of control design discussed in the next chapter requires that 
the representation facilitate efficient manipulation in synthesizing new dynamical 
behaviors. We have chosen an equivalence-class based qualitative representation 
motivated by the above requirements on the form and use of the representation. 

Our goal is to develop a class of intelligent and autonomous programs that 
understand the behaviors of complex [systems through their phase-space represen- 
tations, synthesize control commands, and affect the physical processes. For ex- 
ample, a controller that controls the locomotion of an autonomous walking robot 
would monitor the state of the system, analyze the motion, and command the 
motor to achieve a particular walking behavior. When these systems operate in 
nonlinear regimes and are of high order, their complexities often defy human anal- 
ysis. In order for the robot to autonomously execute control tasks with superior 
performance, it is important that the| robot understands the consequences of its 
control action and synthesizes appropriate control commands to affect the state 
of the system. We are particularly interested in automating the control analysis 

systems that do not lend themselves easily 



and synthesis for a class of nonlinear 



to traditional analysis and design techniques 
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3.2 Qualitative Phase-Space Structures 



We are interested in representing th 
for control analysis and design. One 
space of a dynamical system is in tkms 
stability regions, trajectory flows exhtb 
spatial arrangement of these geometric 
qualitative phase-space structure for 
definition. 



qualitative behaviors of dynamical systems 

useful qualitative representation of the phase 

of equilibrium points and limit cycles, 

iting the same qualitative features, and the 

objects. We call this representation the 

a dynamical system and give the following 



Definition 3.1 Qualitative phases-space structure; The qualitative ph 
space structure within a phase-space 



ase- 



region of interest, for a structurally stable 
dynamical system addressed in this thesis 1 , consists of 

(a) the number, positions, and stability types of equilibrium points and limit cy- 
cles 



(b) the geometric structures of stability regions (basins of attraction) associated 
with the attr actors 

(c) the equivalence classes of trajectory flows 

(d) the spatial arrangement of the equilibrium points, limit cycles, stability re- 
gions, and trajectory flows. 

3.2.1 Equilibria, limit cycjles, and stability regions 

We review some of the basic concepts in dynamical systems theory in order to 
describe the qualitative phase-space structures. Let us consider a single pendulum 



perturbed from its downward resting 
with a smaller and smaller amplitude, 
due to friction. We call such a resting 
an initial state of the pendulum is in 
point if the pendulum starting from 



position. It swings around its vertical axis 
and eventually settles to the resting position 
state a stable equilibrium point. We say that 
the stability region of the stable equilibrium 



that state eventually settles to the stable 
equilibrium point. In this example, all initial states that the pendulum starts from 
except for the vertically upward one a|re in the stability region. 



lr The scope of the class of systems is disclissed in Section 3.5.1 
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Figure 3.1: Equilibrium points: (fr) attractor, (b) repellor, and (c) saddle. 



In general, the equilibrium points of a dynamical system x' = f(x,u), where u 
is a parameter, are the zeros of the vector field f(x,u) : R n — ► R n . Structurally 
stable systems [18] can have equilibrium points of three types: attractor, repellor, 
and saddle, whose local behaviors in phase spaces are shown in Figure 3.1. An 
attractor is an equilibrium point that all nearby trajectories approach in forward 
time. In the pendulum example, the downward resting state is an attractor. We 
call the attractor an asymptotically stable equilibrium point. A repellor is the one 
that repels all nearby trajectories. One can think of it as an attractor in reverse 
time. Trajectories approach a saddle in some directions and leave it in the other 
directions. The vertically upward st^tte of the pendulum is a saddle. Start the 
pendulum resting at this position. A slight perturbation would move it from the 
position. On the other hand, the penjdulum with just the right amount of energy 
could swing to the upward position 4-nd rest there, although it will be extremely 
rare and will take an infinite amount of time to do so. We call the repellors and 
saddles unstable equilibrium points. The study of equilibria plays a central role 
in dynamical systems theory. The equilibrium points, together with other phase- 
space entities, partition a phase spac^ into qualitatively different subregions. The 
attractors are the only kind of equilibria that can be observed physically; they 
characterize the asymptotic behavior^ of dynamical systems. The other classes of 
steady-state behaviors addressed in t lis thesis are limit cycles. 

The important concept of stability is associated with the stability regions. The 
collection of trajectories approaching an equilibrium point is called the stable tra- 
jectories (or stable manifold) of the point; and the collection of trajectories leaving 
an equilibrium point is called the unstable trajectories (or unstable manifold) of the 
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Figure 3.2: A stability region (basin (j)f attraction) for an attractor in a 2nd-order 
system. 



point. A saddle has stable trajectories along some directions and unstable trajecto- 
ries along the other directions. The ujiion of the stable trajectories of an attractor 
is its stability region, often called thje basin of attraction for the attractor. The 
stability region of an attractor is opejn. Every trajectory starting in the stability 
region will be attracted to the attrac|tor, by definition. Attractors cannot be on 
the boundary of the region. The region may be either bounded or unbounded; it is 
unbounded if the boundary contains no repellors. An example of a stability region 
is shown in Figure 3.2. 



3.2.2 Trajectory flows 

Most nonlinear systems do not have; solutions that can be expressed in closed 
form. Analysis of such systems often resorts to computational exploration of their 
trajectories over a period of time. Since a phase space is a union of an infinite 
number of trajectories, an exhaustive search for plausible control paths in the 
space of all the trajectories is computationally too expensive and numerically too 
sensitive to uncertainties. To make it possible to find good control paths, we bundle 
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Figure 3.3: I A trajectory flow. 



ries. Figure 3.3 shows, for example, a 
of all trajectories ending at the same 



together trajectories with similar qualitative features into a trajectory flow. 

The trajectory flows describe the collective shape and direction of the trajecto- 

trajectory flow in the phase space consisting 
destination marked by the symbol +. 

With the trajectory flow bundle representation of the phase space, the explo- 
ration is constrained in a much smaller space of equivalence classes of trajectories. 
The expensive fine- grain search in phase space, possibly of high dimensions, is 
avoided. Furthermore, the flow bundle is a natural representation with respect 
to the robustness to the effect of noise and uncertainties, compared with individ- 
ual trajectories. The effect of noise on an individual trajectory is modeled by 
thickening the trajectory. 

3.3 Automated Qualitative Analysis of Phase- Space Struc- 
tures 



MAPS understands qualitatively different 
metric shape information about these 



regions and extracts and represents geo- 
regions. Given a dynamical system specified 
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as a system of governing equations, 
scription of the phase-space structure 
description will be used as input to 
example, for synthesizing control laws 



MAPS generates a high-level symbolic de- 
as the result of the analysis. The high-level 
other programs for further computations, for 
in the next chapter. 



3.3.1 Theoretical characterization of stability regions 

An essential part of the stability anatysis for a dynamical system is the determina- 
tion of the boundaries of stability regions — the stability boundaries. Our algorithm 
for determining stability boundaries is based on a crucial theoretical result char- 
acterizing the stability boundaries of a fairly large class of dynamical systems 2 . 
Under certain weak conditions to be discussed in Section 3.5.1, the result of Chi- 
ang et al. [10] shows that the stability boundary for an attractor consists of the 
stable trajectories (stable manifold) of equilibrium points and limit cycles whose 
unstable trajectories (unstable manifold) approach the attractor. This allows us to 
numerically determine a collection of trajectory points on the stability boundary 
through calculations of the stable and unstable trajectories 3 . In planar systems, 
the boundaries consist of curve segnjients obtained from trajectories. In higher 
dimensions, the boundary surfaces are swept out by a set of trajectories on the 
boundary. The following theorem establishes the theoretical basis for the stability 
analysis of the algorithm. 

Theorem 3.1 Characterization off stability boundaries for attractors [Chi- 
ang et al. [10]]: For nonlinear autonomous dynamical systems satisfying certain 
conditions, the stability boundary for \an attractor is the union of stable manifolds 
of equilibrium points and limit cycles], whose unstable manifolds approach the at- 
tractor. 



The three conditions that the systems have to satisfy will be discussed later 
in the thesis. The equilibria that can be on the stability boundaries are saddles 
and repellors. A boundary that separates stability regions for different attracting 
sets is called a separatrix. A degenerate boundary is one that does not separate 
stability regions. 



2 Our algorithm does not use Liapunov functions to construct stability regions. The class of 
Liapunov function based methods will be discussed in Section 6.2. 

3 The trajectory points can be obtained, for example, from numerical integration algorithms. 
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Parker and Chua have generalized 
aries for attracting limit cycles 



Figure 3.4(a) shows the stability boundary for the attractor labeled by the 
rightmost +. One of the unstable trajectories of the saddle, denoted by 0, ap- 
proaches the attractor as t — > oo. Therefore, the stable trajectories labeled "b" 
of the saddle form the stability boundary for the attractor. Figure 3.4(b) shows a 
stability region with degenerate boundary labeled "b". 

Theorem 3.1 to characterize stability bound- 
The corresponding result remains much the 
same as Theorem 3.1 except that the occurrences of "attractor" in Theorem 3.1 are 
substituted with the words "attracting limit cycle". Figure 3.5 schematically illus- 
trates the stability boundary for an attracting limit cycle. The boundary is formed 
by two repelling limit cycles. Parker and Chua's generalization is known to be true 
for many important examples. However, it is largely an intuitive generalization; 

Although our techniques for extracting ge- 
attractors are readily applicable to those of 



no rigorous proof has been given yet. 
ometries of stability regions for point 
attracting limit cycles, we shall not discuss them in detail here 



3.3.2 Extracting and representing shapes of stability re- 
gions 

Since a phase space of a nonlinear dynamical system often consists of qualita- 
tively distinct points and regions, the) phase-space "shape" of the system refers to 
the geometric information about the structures and spatial arrangements of these 
points and regions. A qualitative analysis of the dynamical system determines the 
"shape" of the phase space. The geometric information about these regions, for 
example, is extremely useful in analy2;ing stability properties of control designs for 
complex dynamical systems, such as electric power systems and mechanical control 
systems. 



Definition 3.2 The "shape" of a 

namical system refers to the geometric 
points and regions of the phase-space 



We need to extract the geometric 
the numerical results about the stali 
niously so as to facilitate further 
is to be used to estimate the volumes 



dynamical system; The "shape" of a dy- 
characterization of the qualitatively distinct 
representation of the system. 



information about the stability regions from 

lity boundaries, and represent it parsimo- 

conjiputations. For example, the representation 

of the stability regions, to reason about the 
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(a) 




Figure 3.4: Stability boundaries for aiji attractor: (a) non-degenerate boundary 
separatrix; (b) degenerate boundary. 
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attracting, 
limit cycle 




stability boundary 



Figure 3.5: The stability boundary for an attracting limit cycle. 
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spatial relations with the stability b( 
of the regions, to extract geometric 



oun< 



daries, to compute topological properties 
information about trajectory flows, etc. 



Generation of discrete trajectory points 

The boundary of a stability region is numerically approximated by a collection 
of trajectory points. A one-dimensional boundary curve is approximated by a 
sequence of points on boundary trajectories, a two-dimensional boundary surface 
is approximated by a collection of points on boundary trajectories that sweep out 
the surface, and so on. 

The discrete trajectory points are generated with numerical integration method. 
The method of integrating trajectories from initial points is well-know. However, 
it is difficult to generate a set of trajectory points that evenly populate a surface. 
The difficulty arises from deciding wjiere the initial points come from and when 
the integration stops. 

Consider for example a two-dimeijsional stability-boundary surface. Since the 



stability boundary is formed by stable 



manifolds of saddles, the trajectories sweep- 



ing out the surface can be generated by integrating the trajectories backwards from 
initial points near the saddles in the directions of stable eigenvectors. This inte- 
gration method also generalizes to higher-dimensional boundary surfaces. 

In order to maintain relatively even spacings between surface trajectories, the 
program must select a subset from all' the possible initial points and terminate an 
integration of a trajectory when the trajectory gets too close to other trajectories. 
We use the following criteria for this; purpose: (1) the distribution of trajectory 
points, (2) inter-trajectory distances, and (3) local curvature of the surface. 

Let's take the surface of a sphere as a stability boundary, and assume a saddle 
is at the north pole and a source at the south pole. As we integrate backwards from 
the saddle along a set of surface trajectories, these backward trajectories diverge as 
they get near the equator. If the distance between two adjacent trajectories is too 
large, an additional trajectory is needed to cover the gap. The program backtracks 
in this case: it picks an initial point ne&r the north pole, i.e., the saddle, to generate 
this new trajectory. The local curvatjire of the surface also helps determine if an 
additional trajectory is necessary. As the backward trajectories begin to converge 
at the south pole, some of the trajectories can be terminated for the purpose of 
approximation. 



34 



structures on the points, a minimal 
preserving topological structure — is a 



A simplicial representation 

Given a set of trajectory points on the stability boundary, a structure on the points 
is needed to make explicit the proximity relationships among these points on the 
boundary. We define the structure on the trajectory points as a graph: the vertices 
are the given trajectory points; the edges join those points that are related to each 
other, where the notion of relation is (defined by some metric. Among the different 

representation — one with fewest edges and 
polyhedral structure having those boundary 
points as vertices. Furthermore, the resulting polyhedron is contained in the convex 
hull of the boundary points. This polyhedron is not unique. We choose the one 
computed from the Delaunay triangulation to be discussed shortly. 

A polyhedral structure is a consequence of the so-called simplicial represen- 
tation. The elements of a simplicial representation consist of simplices [33]. An 
n-dimensional simplex, or an n-simpjex for short, determined by n + 1 geomet- 
rically independent points 4 is the convex hull of the points. These points act as 
vertices of the n-simplex. A familiar example of a 0-simplex is a point, a 1 -simplex 
is a closed line segment, a 2-simplex is a closed triangle, a 3-simplex is a solid 
tetrahedron, and so on. 

The simplices are basic building blocks from which we construct complicated 
geometric spaces. There are certain ways simplices can be joined to form a new 
object; for instance, two triangles can be glued together at a vertex, or along 
an edge to form a quadrilateral. Given a set of points in a space, a simplicial 
tessellation, or a triangulation, is a collection of simplices that cover the convex hull 
of the given points and have these points as vertices. For example, a triangulation 
on a set of points on a two-dimensionalj surface divides up the surface into triangular 
regions. 

Among possible triangulations over these points, the Delaunay triangulation [40] 
offers certain advantages over the others with respect to the quality of triangles 
measured by the regularity of the triangular mesh. The Delaunay triangulation has 

any simplex does not contain any vertices in 
Delaunay triangulation very attractive: the 



the property that the circumsphere of 
its interior. This property makes the 



?N 



4 A set of points {x , x\, ..., x n ] of R 
and only if the vectors x\ — x , ..., x n — x 
in R 1 form a geometrically independent set 
non-coplanar points in R 3 . 



(N > n) is said to be geometrically independent, if 

are linearly independent. Thus two distinct points 

so do three non-collinear points in R? and four 
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Figure 3.6: A Delaunay triangulation over a set of points in a plane. 



simplices in the triangulation are mos|t equiangular over all triangulations; in two 
dimensions, the minimum interior angles of its triangles is maximized. Figure 3.6 
shows the Delaunay triangulation on a set of points. 

We are concerned with deriving a discrete representation for the geometric 
structure of a phase space that computer programs can efficiently reason about. 
The representation needs to preserve the gross features of shape of the struc- 
ture with minimum complexity. Yip has used a point-set representation for the 
phase-space geometries and topologies in analyzing two-dimensional Hamiltonian 
maps [52]. While the point-set representation has been successfully used to recog- 
nize patterns of the phase portraits, it is not suitable for the purpose of cutting and 
pasting trajectory flows and synthesizing desired phase-space geometry and topol- 
ogy. The simplicial structure, on the other hand, is a much coarser representation 
for describing qualitative features of geometric pieces. The simplicial represen- 
tation is independent of the dimensions, permits recursive decomposition, and is 
characterized by certain algebraic signatures such as fundamental groups [33]. A 
point-set representation or trajectory representation would be too fine-grain for 
manipulating the shape of phase spaces and would incur unnecessary complexity 
in searching for control paths in phase space. 
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Extraction of geometric structures 

The geometric information about a stability region is represented by the polyhe- 
dron tightly stretched over the trajectory points on the stability boundary. Ex- 
traction of this polyhedral approximation proceeds in two steps: computing a 
triangulation of the convex hull containing the polyhedron, and eliminating exte- 
rior triangles. The triangulation is computed with the Delaunay method on the 
set of points which results in a tessellation of the convex hull of these points with 
simplices. The polyhedral approximation is then extracted from the triangulation 
by a sculpture method used in visual information representation [7], followed by a 
more expensive centroid method. Simplices exterior to the polyhedron are elimi- 
nated first by heuristic rules and then by testing their membership in the stability 
region with trajectories emanating from the centroids of the simplices. 

We have implemented the Delaunay triangulation in n dimensions that supports 
incremental insertion and deletion of points, adapted from an algorithm proposed 
for planar triangulation [37]. The time complexity of the algorithm grows exponen- 
tially with the dimensions, so does the number of simplices in the triangulation. 
For N points in n spatial dimensions, the time and space complexities scale as 

n°w [40]. 

Under the condition that the distribution of boundary points is reasonably 
dense and uniform on the boundary, the polyhedron approximating the stability 
region is contained in the triangulation of the convex hull. Furthermore, the union 
of all the circumspheres of interior simplices is a good approximation for the region. 
In order to extract the polyhedron, simplices exterior to the polyhedron have to 
be eliminated. 



For geometric structures with relat 
the exterior simplices are eliminated w 
using the circumsphere heuristics [7] 



tively smooth surfaces, i.e., small curvatures, 
ith a relatively fast and inexpensive method 
We observe that only certain type of sim- 
plices that satisfy the so-called "visibility conditions" can be exterior ones and are 
therefore candidates for elimination. We associate a value V with each candidate. 
The value V is a goodness measure of the polyhedral approximation to the true 
boundary contributed by this simplex. The minimum of all V's of the interior 
simplices is denoted by V min , a measure for the overall goodness of the approxima- 
tion. We focus the search on candidates satisfying the "visibility conditions" and 
eliminate the exterior simplices by deleting the candidate simplex with the least 
V among all the candidates, until the number of points on the boundary of the 
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polyhedral approximation is the same as that of the original set of the boundary 
points and the elimination of this candidate could not further increase V m i n . 

We explore the properties of interior circumspheres versus those of exterior 
ones. When the distribution of boundary points respects the local geometries of 
the boundary they approximate, the envelope of all the circumspheres of interior 
simplices is close to the true boundary. Consider the maximum distance between 
the faces of a simplex that are on the boundary and its circumsphere: the smaller 
the distance is, the better the simplex approximates the boundary. The distance 
for an interior simplex is smaller than that of an exterior one. We use the inverse 
of this distance as the value of V for the simplex. Another possible candidate 
for V is the inverse of the circumsphere radius. Figure 3.7(a) illustrates that for 
smooth shapes the exterior triangles tend to be flat and have larger circumcircles, 
compared with the interior ones. 

We reiterate that the condition oil the distribution of the boundary points has 
to be checked with respect to the shjape of a region, to ensure that the circum- 
sphere heuristic rule for elimination A|vorks. For example, more points are needed 
to approximate the boundary of a region with finger-like narrow parts, as shown 
in Figure 3.7(b). I 

The sculpture method is stated as follows: 

Algorithm 3.1 Sculpture method [7]: 

1. Construct a set G of all simpliced in the triangulation and a set H of candidate 
simplices satisfying "visibility conditions" ; 

2. Order the simplices of H according to their V values; 

3. Consider the simplex s of H having the least V. If the number of vertices of the 
polyhedron formed by the set G is less than that of the original set of points or 
the removal of s increases the V nuri of H , then delete s from G and H] otherwise, 
return G and terminate; 

4. Update neighbors n t - of s in G. Add ni to H . Go to Step 2. 



Visibility conditions: 

• In 2D, simplices with one edge piid two vertices on the boundaries. 
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smooth shape 



(a) 



narrow part 




(b) 



Figure 3.7: Circumcircle properties 
region with smooth shape; (b) a regi 



of exterior triangles vs. interior ones: (a) a 
with narrow-parts. 



ion 



39 



• In 3D, simplices with exactly oite 
boundary, or those with exactly 
boundary. 



3 face, three edges, and three points on the 
two faces, five edges and four points on the 



• In n-dimensions, simplices with one, or two, ..., up to n — 1 proper faces on 
the boundary. ■ 

Each simplex s in the triangulation has a list of neighbors that share one or 
more common faces with s. When s is ideleted from the triangulation, its neighbors 
need to be updated. A simplex satisfying the "visibility conditions" is added to 
H once and is deleted from H at most once in the elimination process. Let the 
number of simplices in the set G be M. When the set H is implemented as a heap, 
each insertion or deletion of H costs al; most 0(log M). The worst-case complexity 
of the algorithm therefore is 0(M log M). 

For systems whose stability regions have interleaved geometric structures or 
holes, the geometric extraction algorithm using the circumsphere heuristic de- 
scribed earlier often terminates before| all the exterior simplices are eliminated. In 
this case, MAPS resorts to a more expensive method, the centroid method, that 
determines the membership of each simplex in the stability region: a simplex is 
classified as being in a stability regioiji of an attractor if the trajectory emanating 
from the centroid of the simplex approaches the attractor in the limit or enters 
another simplex already in the stability region. For example, the buckling column 
discussed in Chapter 2 exhibits a bajided phase space (see Figures 2.1 and 2.4). 
The circumsphere heuristic method terminates before all the boundary trajectory 
points are on the boundary of the polyhedral approximation. The centroid method 
is required to eliminate the remaining triangles. 

The centroid method is implemented by an unraveling algorithm, adapted from 
the cell-to-cell algorithm [22], which links together simplices according to the map- 
ping of their centroids under trajectory flows. 

Algorithm 3.2 Centroid method: 

1. Let a be the attractor under consideration, A be the set of simplices containing a in 
the interior or on the boundaries, G be the set of all simplices in the triangulation 
except for those of A, and C be a set initially empty; 
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2. For a simplex s G G\ let p be the 
the centroid of s enters. Delete s 
simplices 5 from s to A; otherwise, 



next simplex that the trajectory starting at 
from G. If p G A, add s and all the reachable 
add s to C and make a pointer from p to s; 



3. If G is empty, return A and terminate; otherwise, go to Step 2. 

The time complexity of centroid hiethod scales linearly with the number of 
simplices in the triangulation. The method requires integration over a time interval 
only once for each simplex. 

The sculpture method and the centroid method work together to extract the 
geometric shape of the stability regions. The information about the boundary is 
then used in constructing the flow-pipe representation for the phase space. More 
specifically, the boundary information is used in labeling the simplices with flow 
directions at the boundary. Details of this will be discussed next. 



3.3.3 Modeling trajectory flow pipes 

The flow pipes are further geometric and dynamical characterizations of stability 
regions. They form the foundation fori the search algorithm for synthesizing global 
control paths. 

The polyhedral modeling of stability regions provides a geometric structure 
for reasoning about the stability boundaries. However, it does not provide in- 
formation on how the trajectories floiv within the region — the transient behavior 
very important in high-quality control design. The characterization of trajectory 
flows captures the missing information about the underlying vector-field flows; 
such information is essential for efficiently "looking for" good control reference 
trajectories. 

We introduce the technique of floW pipes for describing the direction and the 
shape of the trajectory flows. The flow pipes form a discrete representation for 
phase space and yet preserve the fine geometries of the trajectories. The flow- 
pipe modeling, together with that of equilibria and stability regions, characterizes 
the transient and asymptotic behaviors of the dynamical systems in terms of the 
qualitative phase-space structures. Figure 3.8(a) shows a portion of the phase 
space for the Lienard equation x" -f 0.5x' -fa; 2 = 1, consisting of all trajectories 
ending at the attractor denoted by the symbol -f. A flow pipe groups this collection 



5 A simplex x is reachable from s if there 



is a pointer or a chain of pointers from s to x. 
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of trajectories that exhibit the same Qualitative feature into an equivalence class, 
as shown in Figure 3.8(b). 



models a collection of trajectories exhibit- 
e phase space. It is an equivalence class of 
continuously deformed to another one, de- 
regions, trajectories connecting equilibria, or 



Definition 3.3 flow pipe: A flow pipe 
ing the same qualitative features in th 
trajectory paths, each of which can be 
limited by the boundaries of stability 
boundary surfaces dictated by applications" 

Constructing flow pipes 

Computationally, a flow pipe is constructed by aggregating consistent geometric 



pieces with respect to the way that 
catalogue of flow types on boundaries 



the flow travels from piece to piece. The 
of the geometric pieces is a set of consistency 



constraints derived from the underlying dynamics of the vector field 

Given the simplicial tessellation of the phase space constructed by MAPS, each 
proper face / of a simplex s in the tessellation is classified into one of the three 
types with respect to the direction of the flow (see Figure 3.9), according to the 
following labeling scheme: 

1. unidirectional flow: 

label in: Trajectories flow into s on /; 
label out: Trajectories flow out of 5 on /. 

2. bidirectional flow: 

label in-and-out: Trajectories flJDW into s on some portion of / and flow out 
of ^ on other portion. j 

3. zero flow: 

label tangent: The flow does no|t intersect /. Hence the proper face / is on 
the boundary of the flow. 

To determine the direction of th£ flow on a face of a simplex, the program 
samples N points evenly distributed on the face and computes the components of 



6 Formally, a flow pipe models a homotbpy 
by the boundaries of stability regions, trajectories 



boundaries. If the ends of a flow pipe are 
equivalence class of trajectories. 



equivalence class of trajectory paths delimited 
connecting equilibria, or application-specific 
squeezed to points, the flow pipe is a path homotopy 
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(a) 




Figure 3.8: Grouping trajectories into 
with the same qualitative feature; (b) 



(b) 



flow pipes: (a) a collection of trajectories 
a corresponding flow pipe. 
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in-and-out 
(2) 



tangent 
(3) 



Figure 3.9: Classifying faces 



of simplices with respect to flows. 



flow vectors along the normals of the face at these points. The stability boundaries 
necessarily separate flows and hence are on the flow-pipe boundary. Since the true 
stability boundary is only approximated by polyhedral faces, the above sampling 
method could yield incorrect flow labels on the boundary. The boundary infor- 
mation obtained in the characterization of the stability regions described earlier 
is used in labeling these zero flow boundary faces. In addition, the trajectories 
that connect saddles on the boundary and the attractor also separate trajectory 
flows. The simplicial tessellation is refined with the introduction of trajectory 
points sampled along these connecting trajectories. 

A proper face of a simplex is called monotonic with respect to both the flow 
and the simplex, if the face is labeled either type 1 or type 3 in the above labeling 
scheme. A simplex is monotonic with respect to the flow if all of its proper faces 
are monotonic. Similarly, a polyhedron is monotonic if all of its proper faces are 
monotonic. Given two simplices agreeing on a non-monotonic proper face, the 
common non-monotonic face is canceled in the polyhedron formed by the two 

-monotonic polyhedron continues to grow in 
are monotonic. The monotonic polyhedra 



simplices; see Figure 3.10(a). A non- 
this way until all of its proper faces 



are then aggregated to form flow pipes in such a way as to conserve the flow on 



their proper faces; see Figure 3.10(b} 
into smallest monotonic polyhedra an 
pipes. 



The following algorithm groups simplices 
d clusters the monotonic polyhedra into flow 



Algorithm 3.3 Flow-pipe construction: 
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Figure 3.10: Construction of flow pipes 
polyhedron by canceling common non 
polyhedra to form a flow pipe accordi q 



: (a) grouping simplices to form a monotonic 
monotonic faces; (b) aggregating monotonic 
g to flow directions at the boundaries. 
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1. Label proper faces of each simplex 



2. Cluster simplices into monotonia 
Cluster simplices into equivalence 
proper faces of the simplices: twc 
face with the same flow label: in- 



with respect to the flow. 



polyhedra: 

classes with the equivalence relation on the 
simplices are equal if they agree on a proper 

dnd-out. 



3. Cluster monotonia polyhedra into flow pipes: 
Cluster monotonic polyhedra into equivalence classes with the equivalence relation 
on the proper faces of the polyhedra: two monotonic polyhedra are equal if they 
agree on a proper face with consistent flow directions: in and out, respectively. 

4. Order the monotonic polyhedra i^ each class from Step 3 according to the flow 
direction. Return. 



the 



The resulting flow pipes model 
pipe is formed by the proper faces labeled 



Classifying shapes of flow pipes 



net 



an 



Flow pipes capture dynamical shapes, 
starts from a repellor. Saddles can 
of the pipe, since flow pipes further 
region, each flow pipe either ends at 

Flow pipes exhibit topological 
its relative shape. A flow pipe aroun 
pipes are open. The open pipes can 
two, or more pipes can side by side 
Figure 3.11. We classify the pipes inti 



shapes 



win 



; trajectory flows. The boundary of a flow 
tangent of simplices in the pipe. 



A flow pipe ends at an attractor and/or 
be in a flow pipe except for the boundary 
decompose the stability regions. In a bounded 
attractor or leaves the region. 

Each flow pipe is classified according to 

a stable limit cycle is a closed pipe. Most 

e relatively straight or highly wound. One, 

d into a spiral, called n-winding pipes. See 

the following categories: 



1. open pipes: (a) straight pipes; (b) wound pipes. 

2. closed pipes. 



In the earlier Figure 3.8(b), for exa|mpl 
at its end and forms the stability 
patterns can occur with more pipes or 



e, a single pipe winds into a spiral shape 
region for the attractor. More complicated 
in higher dimensions. 
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straight 



(b) closed pipes 



Figure 3.11: Topolo, 



gical shapes of flow pipes. 



3.4 The MAPS Analysis jAlgorithm 
3.4.1 The algorithm 

We present the following algorithm for analyzing, extracting, and representing 
qualitative features of a dynamical system of any order in the phase space. 

Algorithm 3.4 MAPS analysis: 

(1) Identify qualitative behaviors: 

(a) locate equilibrium points/limit cycles and classify their stability types; 

(b) compute stable and unstable trajectories for each saddle/limit cycle; 

(c) identify those saddles/limit cycles 
an attractor; 

(d) the stability boundary for the attfactor is the union of the stable trajectories 
of those saddles/limit cycles identiied in Step (c); 

(e) check if consistency rules are violated. If yes, look for missing equilibrium 

(a). Otherwise, go to the next step. 



points/limit cycles and go to Step 



(2) Extract geometric structures: 



whose unstable trajectories approach 
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(a) for each attractor, collect stabilit 

(b) tessellate the convex hull of the b 

(c) extract a polyhedral approximation 



boundary points; 
oundary points with a tri angulation; 
to the stability region. 



(3) Construct flow pipe: 

(a) refine the triangulation of the stability 

(b) label simplices with flow direct 

(c) aggregate the labeled simplices kit 



regions; 
ions at boundaries; 
o flow pipes. 



(4) Summarize qualitative behaviors 

(a) compile the phase-space data stri cture 

(b) augment the graph with the geometric 

(c) report the graph as the output. 



The set of consistency rules specify 
and are used in the algorithm to automatically 



1. The Existence Rule: Every stability 
phase space with multiple attractors 



2. The Separation Rule: Separat 
bounded on all ends. 

The first rule states the existence 
with multiple attractors. The second 
multiple stability regions. The separat 
two stability regions. 

The first step of our algorithm is 
Parker and Chua [38] for numerical 
nar systems. We have augmented 
they suggested to automate the locating 
method [41] used in finding equilib 
Chua method uses a grid to set up ini 
and unstable equilibrium points und 
quire that the initial guesses for sad 
seek to automate saddle locating by 
most likely places using partial boundary 
fining the initial guesses for the Newton 



their 



er 



and generate a high-level description: 

from Step 1 into a relational graph; 
: structures from Steps 2 and 3; 



the conditions for the stability boundaries 
locate missing saddles. 



region of an attractor has a boundary in a 
nc0s either form a closed surface or become un- 



of stability boundaries in a phase space 

rule describes the separation property of 

rices are stability boundaries that separate 



based on a numerical method proposed by 

y determining stability boundaries of pla- 

ir method with the set of consistency rules 

of saddles. Since the Newton-Raphson 

points requires an initial guess, the Parker- 

tial guesses and is able to find all the stable 

normal circumstances. However, they re- 

es be provided manually by the user. We 

ing the search for missing saddles on the 

information already obtained, or by re- 

Raphson method. We want to emphasize 



dl 



focus 
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INPUT 



system 



Locate equilibriuin 
and classify 



local 



Compute saddle 
determine s 



trajectories and 
tability regions 



Extract 

of stability 



Summarize 
phase 



equation 



points/limit cycles 
stability types 




geonjietric structures 
regions 



Construct flow pipes 



space 



and represent 
structures 



HIGH LEVEL DESCRIPTION 

relational graph 



Figure 3.12: The flow chart of MAPS. 
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ace 



an ax 



there 



sequence 



saddl* 



that our algorithm is valid for higher 
symbolic description of the phase-sp 
signed for numerically analyzing pi 
also finds degenerate boundaries if 

In the second step, the algorithm 
ity regions. It uses a triangulation 
system and computes polyhedral 
represents a trajectory with a 
4th-order Runge-Kutta method [41] 
boundary points that approximate a 

The third step refines the trian 
flow pipes, using the information 
trajectories of saddles connect the 
jectories the connecting trajectories 
the insertion of trajectory points from 
for the flow pipes clusters simplices i 
tion 3.3.3. 

The last step represents the data 
topologies with a relational graph 
an example in the following section. 

The flow chart of MAPS is show 
system of governing equations for a 
a set of measured states from e 
be analyzed and clustered to form p 
be reconstructed from a measured 
technique called delay coordinates 



dimensional systems as well and generates a 
; structure. The Parker- Chua method is de- 
systems only. We note that our algorithm 
are any. 
extracts the geometric structures of stabil- 
method to tessellate the phase space of the 
approximations to the regions. MAPS internally 
; of integration points. It uses an adaptive 
to integrate trajectories and to generate the 
tability boundary, 
gr.lation from the second step and constructs 
about the stability boundaries. The unstable 
es with the attractors. We call these tra- 
The triangulation from Step 2 is refined by 
the connecting trajectories. The algorithm 
Into equivalence classes as described in Sec- 



He 



3.4.2 The main illustration 

We illustrate how MAPS computes 
with an example. Consider a 2nd- 



J x' = —3x -f 4x 2 — xy/2 — 

\ y' = -2.1j[ 
where u is a parameter. For the parameter 
the region — 1.0 < x < 4.0 and —1.0 



structure of the phase-space geometries and 
; details of each step will be illustrated with 



in Figure 3.12. The input to MAPS is a 

dynamical system. We could also start with 

experiments. These numerical states could then 

base portraits. The phase space could also 

series of just one state variable, using a 

introduced in [36, 48]. 



time 



the high-level description of a dynamical system 
order nonlinear system 



(3.1) 



+ xy + u 

value u = 0.2, the vector field within 
< y < 4.0 is shown in Figure 3.13. MAPS 
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Figure 3.13: The main illustration: the vector field of a 2nd-order nonlinear system. 



ginal 



analyzes the qualitative behavior of 
actually performs the analysis in a 
trajectories that enter the bounding 
implementation magnifies the ori 
side.) The equilibrium points of the 
on f(x,u) — the Newton-Raphson me 
within the region and classifies their 
Jacobians at the equilibrium points: 
and two saddles at (1.05,0.19) and (3 
that the attractors are represented 
sented with the symbol ©.). The stab 
then computed by integrating the syst 
in the directions of the stable and 
respectively. 

Since one of the unstable 
(2.0,2.0), the stability boundary of 
of both saddles. Similarly, the stabili 



wi 



the 



lity 



th 



e system within this region. The program 
lightly larger bounding box to account for 
box again after a short exit. (The current 

bounding box by a factor of 1.5 on each 

system are found by a zero-finding method 

■thod. MAPS locates four equilibrium points 

stabilities by inspecting the eigenvalues of 
two attractors at (0.0,0.0952) and (2.0,2.0) 
f)5, -0.21), all shown in Figure 3.14(a) (Note 
th the symbol + and the saddles are repre- 
e and unstable trajectories of the saddle are 
em from a small neighborhood of the saddle 
e eigenvectors backwards and forwards, 



unstabl 



trajectories of each saddle goes to the attractor at 

attractor consists of the stable trajectories 

boundary of the attractor at (0.0, 0.0952) 
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the 



of 



different 



consists of the stable trajectories of 
trajectories goes to the attractor. H 
trajectories that leave the bounding 
thought of as the stable trajectories 
ble trajectories of the saddle at (3.05 
attractor at infinity, for one of the 
bounding box. Consistency rules are 
MAPS finds three qualitatively 
tors and internally represents the ph 
attractors are connected with each o 
boundaries (Figure 3.14(b)). 

The second step extracts a po 
preserving the gross features of the 
region of the attractor at (2.0, 2.0). 
imated by a collection of trajectory 
uniform and dense on the boundary; 

A Delaunay triangulation is 
convex hull of the points is tessellated 
boundary points are distributed reas 
the polyhedral shape of the stability 
convex hull. MAPS successfully eli 
with the heuristic method, without 
region is shown in Figure 3.16(e). 

MAPS then computes the flow- 
lation of the stability region obtained 
from connecting trajectories. MAPS 
tion with flow direction on its faces 
classes. This results in two flow 
second step (Figure 3.16(f)). 

MAPS compiles the data struct 
pipes from earlier steps into a 
findings to the user in a symbolic 



unst 



The 
points 



saddle at (1.05,0.19), one of whose unstable 

owever, within the region of interest there are 

box. These trajectories can be conveniently 

an attractor at infinity. Therefore, the sta- 

, —0.21) form the stability boundary for the 

table trajectories of the saddle leaves the 

^hecked and satisfied. At the end of this step, 

regions associated with the three attrac- 

ase-space structure in a data structure: the 

ther via saddles and associated with stability 

lyhedral approximation to each stability region 

shape of the region. Consider the stability 

stability boundary is numerically approx- 

s (70 in total for this example) relatively 

see Figure 3.15(c). 

performed on this set of points. As the result, the 

with triangles; see Figure 3.15(d). Since the 

onably densely and uniformly on the boundary, 

region is contained in the triangulation of the 

liminates the triangles exterior to the polyhedron 

reporting to the centroid method. The resulting 

pirie modeling of the phase space. The triangu- 

above is further refined with additional points 

labels each triangle in the refined triangula- 

and clusters these triangles into equivalence 

pipes that form the stability region from the 



ure 



for equilibria, stability regions, and flow 
relatibnal graph; see Figure 3.17. It reports its 
summary in Figure 3.18. 
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Figure 3.14: The main illustration — 
boundary and connecting trajectories 



(a) 




(b) 



MAPS output: (a) equilibrium points; (b) 
(to be continued) 
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Figure 3.15: (cont'd) The main illustrat 
stability boundary for one of the attralctors 



(c) 




(d) 



ion — MAPS output: (c) points on the 
; (d) triangulation of the convex hull. 
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Figure 3.16: (cont'd) The main illustration 
proximation to the stability region comput 
points; (f ) two flow pipes computed from 
form the stability region. 



(e) 




(f) 



— MAPS output: (e) polyhedral ap- 
ed from the triangulation on boundary 
the refined triangulation. The flow pipes 
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stability region 3 

flow pipe 4 

boundary trajs 1,2 

connecting traj 5 




stability region 2 
flow pipes 1,2 
boundary trajs 1,2,3,4 
connecting trajs 6,7 



stability region 1 

flow pipe 3 

boundary trajs 3,4 

connecting traj 8 



Figure 3.17: The main illustration: the relational graph constructed by MAPS. 
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<equilibrium-points : 
equilibrium 1. (saddle at (3.05 -.2i)) 
equilibrium 2. (attractor at (2. 2.}) 
equilibrium 3. (saddle at (1.05 .19]) 
equilibrium 4. (attractor at (0. . 0952))> 



<trajectories : 

<boundary-traj ectories : 
trajectory 1. (from *infinity* t 
trajectory 2. (from *infinity* t„ 
trajectory 3. (from *infinity* to 
trajectory 4. (from *infinity* to 

<connecting-traj ectories : 
trajectory 5. (from (1.05 .19) to 
trajectory 6. (from (1.05 .19) to 
trajectory 7. (from (3.05 -.21) t 
trajectory 8. (from (3.05 -.21) t 



o ( 



<stability-regions : 
stability-region 1. 

attractor at *infinity* 

stability-boundary: (trajectory 4. 

connecting-traj ectories : (traj ectdry 
stability-region 2. 

attractor at (2. 2.) 

stability-boundary: (trajectory 4 

connecting-traj ectories : (traj ectojry 
stability-region 3. 

attractor at (0. .0952) 

stability-boundary: (trajectory 2 

connecting-traj ectories : (traj ectotry 



<f low-pipes: 
flow-pipe 1 . 

boundary : 
flow-pipe 2. 

boundary : 
flow-pipe 3. 

boundary : 
flow-pipe 4. 

boundary : 



(from *infinity* to (2 
(trajectory 3. traject 

(from *infinity* to (2. 
(trajectory 4. trajectory 

(from *infinity* to * 
(trajectory 4. trajecto 

(from *infinity* to (0 
(trajectory 2. trajectorly 



;ory 



=inf 



ry 



1.05 .19)) 
1.05 .19)) 
3.05 -.21)) 
3.05 -.21))> 

0. .0952)) 

2. 2.)) 

(2. 2.)) 

♦ infinity* )» 



trajectory 3. ) 
8) 



trajectory 3. trajectory 2. trajectory 1.) 
7 trajectory 6) 



trajectory 1 . ) 
5)> 



2.)) 

1. trajectory 7. trajectory 6.) 
2.)) 

2. trajectory 7. trajectory 6.) 
Unity*) 

3. trajectory 8. ) 
0952)) 
1. trajectory 5.)> 
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Figure 3.18: The main illustration — MAPS output: a symbolic summary. 



3.4.3 Implementation details 



object 



y 



MAPS locates an attractor at (1.06, 
within the region —5.0 < x < 5.0, - 
determines that the stable trajectories 
for the attractor. The stability boundary 
proximated by a set of relatively evenly 
the phase space with tetrahedra and 
stability region of the attractor (see 
general geometric primitives — the n 
in n dimensions. 

We have also run MAPS on a syst 
domain of our method (see discussion 
2nd-order example (3.1) in Section 3 



MAPS is implemented in the Scheme 
is a vector of the coordinates. An equil 
ues and eigenvectors of Jacobian at t 
specifies the type of the trajectory, s 
gration points as a Scheme stream 
or unstable-man, standing for stable 
typed stable-man is integrated backward 
is integrated forwards. A stability 
the boundary and connecting trajectories 
flow pipe is a list of its polyhedral < 
and end faces. Table 3.1 summarizes 
The examples for the stability region 
limitations. 

3.4.4 More examples 

Our algorithm also applies to higher-dimensional systems. Consider the following 
3rd-order nonlinear system 



Language. Internally, a point in phase space 
ibrium point is a list of its position, eigenval- 
le point, and its stability type. A trajectory 
art and end points, and a sequence of inte- 
The type could be either stable-man 
manifold or unstable manifold; a trajectory 
s, and a trajectory typed unstable-man 
is represented as a list of the attractor, 
, and the polyhedral approximation. A 
approximation, the boundary, and the start 
the internal representation of these objects, 
and the flow pipe are omitted due to space 



x -y 
- y 2 - Q.Sz 2 



0.0,1.12) and a saddle at (-1.06,0.0,1.12) 

5.0 < y < 5.0, and -5.0 < z < 5.0, and 

of the saddle form the stability boundary 

is a two-dimensional surface and is ap- 

spaced trajectories. MAPS then tessellates 

extracts a polyhedral approximation to the 

; 3.19). As described earlier, MAPS uses 

^implices — to approximate stability regions 



em that is not generic and falls outside the 
in Section 3.5.1). The system is the earlier 
4.2 with parameter value u = 0. MAPS 
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Object 

Representation 

Example 



Object 

Representation 

Example 



Object 

Representation 

Example 



Object 
Representation 

Example 



Object 

Representation 

Example 



a point 

a vector of coordinates 

#(2. 2.) 



an equilibrium 

a list: (point, ei 

(#(2. 2.) (-. 

(#(1. 0) #(. 



gerval 



ues/eigenvectors, stability type) 
0,5+1. 41i -.05-1.41i) 
1.)) stable) 



C25 



a trajectory 
a list: (type, start, 
(stable-man *i 
(#(3.05 -.21) 



a stability region 
a list: (attractor, 



polyhedron) 



a flow pipe 

a list: (polyhedron. 



end, integration points) 
infinity* #(3.05 -.21) 
# [promise 26])) 



boundary trajectories, connecting trajectories, 



boundary, start faces, end faces) 



Table 3.1: The internal representation of phase-space data objects. 
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Figure 3.19: The analysis of a 3rd-ordir 
ity boundary and connecting trajectories 
approximation in x-z plane. 



(b) 



nonlinear system: (a) projection of stabil- 
in x-z plane; (b) projection of polyhedral 
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locates four equilibrium points with 
4.0 and computes the stable and 
with the following partial description 



in 



uns 



<equilibrium-points : 
equilibrium 1. (saddle at (3. 0.)) 
equilibrium 2. (attractor at (2.1 1 
equilibrium 3. (saddle at (1. 0.)) 
equilibrium 4. (attractor at (0. 0. 



• 98)) 
))> 



<trajectories : 

trajectory 1. (from *infinity* to fl. 0.)) 

trajectory 2. (from *infinity* to ;1. 0.)) 

trajectory 3. (from (1. 0.) to (3. 0.)) 

trajectory 4. (from *infinity* to (3. 0.)) 

trajectory 5. (from (1. 0.) to (0. 0.)) 

trajectory 6. (from (3. 0.) to (2.:. 1.98)) 

trajectory 7. (from (3. 0.) to *infinity*)> 



<saddle-connection: 

trajectory 3. (from (1. 0.) to (3. 



MAPS discovers a saddle connection 
graph of the phase-space structure: 
trajectory of the saddle at (1.0,0.0) 
(3.0,0.0). The saddle connection is 
concludes that the system is not 
for the planar system here — and ab 
phase-space structure. Since saddle 
structural instability, we believe that 
the phase-space structures of chaotic 



generic 



>an 



3.4.5 Hierarchical 

space information 



We have described and demonstrated 
system through successive computations 
equation representation. MAPS generates 
ical system at the end of the analysis 



the region -1.0 < x < 4.0 and -1.0 < y < 
stable trajectories of saddles. It terminates 



o.))> 



th 



labeled 



in the course of determining the relational 

e trajectory that serves both as an unstable 

and as a stable trajectory of the saddle at 

as trajectory 3 in Figure 3.20. MAPS 

which also implies structural instability 

dons any further efforts to characterize the 

connections are often precursors to chaos or 

they are important in partially characterizing 

or structurally unstable systems. 



extraction and representation of phase- 



our algorithm for analyzing a dynamical 

on the system, starting from its system 

a high-level description of the dynam- 

To bridge the large semantic gap between 
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Figure 3.20: A 2nd-order 



system with a saddle connection. 



the deduced symbolic description and 
put, MAPS has employed multiple in 
system, shown in Figure 3.21. 

MAPS extracts the information i 
each intermediate representation 
properties, such as spatial relations, 
explicit and thus can be accessed an 
level of the representation. In the 
description of equilibrium points, 
computes the polyhedral approximations 
Finally, it pieces together the local 
global picture of the phase-space 
interactions of equilibrium points, 

The internal representation of the 
captures the qualitative aspects of 
graph, nodes are attractors and arcs 
regions. Each node has information 
ated stability region, trajectory flows 



the system equation representation of the in- 
termediate representations for the dynamical 



incrementally, applying a set of operations to 

each level of the representation, implicit 

of the system at different scales are made 

I manipulated by the operators at the next 

of analysis, MAPS first generates a local 

cycles, and their eigenstructures. It then 

to stability regions and trajectory flows. 

information about each limit set to form a 

the relational graph describing the 

cycles, stability regions, and flow pipes. 

hase-space description as a relational graph 

e phase-space structure. In the relational 

denote the relations between their stability 

about the attractor it represents, the associ- 

and their polyhedral approximations, and 



At 



orcer 
limit 



stricture: 



limit 



th 
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INPUT 

system equation representation 



Local dynamical 
(equilibrium points/limi 



(geometric modeling of 



Global 
(interaction of 



Geometric description 

stability regions/flow pipes) 

relational description 
stability regions/flow pipes) 



Symboli 



Figure 3.21: A multi-layered representation for a dynamical system used in MAPS. 



the boundary trajectories and boundalry 
equilibrium point has information about 
and stability type. 



3.5 Discussion 

MAPS analyzes qualitative behaviors 
edge about stability and trajectory 
tracts the geometric information from 
a qualitative phase-space structure, 
ical systems to which our method 
role of MAPS in assisting control 
systems. 



appl 



3.5.1 Scope of the analysis 

We have stated that the theoretical 
boundaries is the result of Chiang et 



♦ 



description 
t cycles, local stabilities) 



♦ 



c description 



equilibrium points and limit cycles. Each 
its position, eigenvalues and eigenvectors, 



of nonlinear dynamical systems using knowl- 
from dynamical systems theory. It ex- 
the numerical results and represents it with 
this section, we discuss the class of dynam- 
ics, further extensions to MAPS, and the 
in visualizing behaviors of dynamical 



flows 



In 



engineers 



basis of our algorithm for locating stability 
al. [10]. Chiang et al. give a topological and 
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dynamical characterization of the stub 
dynamical systems and show that the 
by the stable trajectories of equilibrium 
This class of dynamical systems is defined 



1. Hyperbolicity: All equilibrium pjoints and limit cycles on stability boundaries 
are hyperbolic 7 . 



2. Transversality: The stable and 
and limit cycles on stability b 



unstable trajectories of equilibrium points 
oundaries are transversal to each other. 



3. Finite limits of boundary trajectories 
boundaries approaches an equil 



Conditions (1) and (2) are generic 
cal systems. The example in Section 
condition (2). Structurally stable 
(2) (see [19, 46] for a definition of 
structural stability of dynamical syst 
gineering applications are structurally 
turally stable systems for the sake of 
absence of quasi-periodic orbits. Con 
systems. Our method applies to a 
terns as defined by the three condit 
dynamical system can be converted i 
be analyzed by MAPS. 

In its current implementation, 
attractors, in which stability regions 
structures. It also assumes that the 
system is finite. 



syst 



generic 



ems 



fairly 



ions 



into 



field 



th 



7 Consider the Jacobian of the vector 
the Jacobian have non-zero real parts, then 
cycle is hyperbolic if none of its characteristic 
lies on the unit circle in the complex plane 

8 This is not to say that structurally uns 
systems are important in investigating ph 



stabl 



ility boundaries for a class of autonomous 
stability boundary of an attractor is formed 
points and limit cycles on the boundary, 
by three conditions: 



s: Every trajectory on the stability 
brium point or a limit cycle as t — » oo. 



• properties satisfied by almost all dynami- 
3.4.4 that has a saddle connection violates 
ems satisfy the generic conditions (1) and 
properties and [18] for discussions on 
). Since many systems encountered in en- 
stable, we further restrict ourselves to struc- 
simplicity 8 . Structural stability implies the 
dition (3) excludes some structurally stable 
arge class of autonomous dynamical sys- 
. We note that a periodic non-autonomous 
an autonomous system [18] and can also 



MAPS does not handle systems having chaotic 
and/or stability boundaries exhibit fractal 
number of equilibria and limit cycles of a 



at an equilibrium point. If all the eigenvalues of 

e equilibrium point is hyperbolic. Similarly, a limit 

multipliers, the eigenvalues of its Poincare map, 



e systems are not interesting. In fact, Hamiltonian 
in which energy dissipation is negligible. 



enomena 
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>ensive 



ilat 



3.5.2 Extensions 

The theoretical basis of our approach 
ality, as does the geometric extraction 
approximation is a natural, economical 
We have run MAPS on systems in 
the computation can be very exp 
complexity of the Delaunay triangu 
dimensions. The number of points 
a function of the quality of the 
reflected in managing the large data 
tries needs to be addressed at greater 
a reduction of the complexity is possi 
strained to a submanifold of a high 
recognize such structure of the phase 
subspace. 

In dimensions three or higher, a 
jectories sweeping out the surface 
trajectories of the boundary to 
currently uses spacings between adj 
ture to measure the quality of the 
other kind of metric information abdut 
approximation. Much remains to be 
plicated curved hypersurfaces and to 

MAPS extracts the geometric 
bility boundaries and additional 
plicated stability regions, such as th 
this approach. More work needs to b 
shapes with various kinds of approaches 

We have shown that MAPS detected 
ble system. A considerable amount of 
robust methods for detecting saddle 
nize chaotic attractors. In the case wh 
the fine structure of the boundary 
surface, it is still possible to envelop 
certain thickness. Another possible 



dimensions 



lose 



holds regardless of phase-space dimension- 
method of our algorithm. The polyhedral 
representation for the boundary surface, 
up to three. In higher dimensions, 
We have noted in Section 3.3.2 that the 
ion grows exponentially with the number of 
necessary to approximate the boundary is also 
approximation. The computational complexity 
structure for representing phase-space geome- 
length in future research. One case in which 
ble is when the interesting dynamics is con- 
dimensional phase portrait. A program could 
space and restrict the computation to that 

stability boundary surface consists of tra- 
Obtaining a set of relatively evenly spaced 
approximate the surface is challenging. MAPS 
acent boundary trajectories and local curva- 
appjroximation. Combining these measures with 
the boundary surface could improve the 
explored in order to best approximate corn- 
reflect the underlying dynamics, 
structure using the information about the sta- 
infoifmation about the trajectories. Other com- 
; containing holes, can also be tackled with 
done to catalogue the extraction of region 



a saddle connection in a structurally unsta- 

work needs to be directed toward exploring 

connections and extending MAPS to recog- 

'311 the stability boundary is fractal, although 

would be difficult to characterize with a crisp 

e fractal structure with a boundary layer of 

extension is to augment the current program 
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with a bifurcation analysis, similar to 



3.5.3 The use of MAPS in visualization 

can assist engineers in designing controllers 



The qualitative phase-space descript 
for complex systems. 

We envision that a control eng 
of phase-space structures when varyin 



ion 



meer 



1. the birth, disappearance, and njovement of equilibrium points and limit cy- 
cles; 

2. the enlargement and shrinking of stability regions; 

3. other phase-space qualitative changes. 



rap hi 



:rai 



past* 



A program has been constructed for 
metric structures: it renders the surfaces 
lighting parameters. With such a g 
alize the geometries of a phase port 
corresponding changes, and cut-and- 
different parameter values to form a 

In the next chapter, we automat^ 
Navigator that autonomously synth 
More importantly, the method will b 
tern whose phase space is high dimen 
the desired control properties are impcfssibl 
niques. The topological and dynamical 
and trajectory flows forms the basis 



3.6 Summary of the Chapter 



We have developed a qualitative meth 
structures of nonlinear dynamical syst 
strate the method. MAPS "looks" 
regions — the stability regions, models 
extracts and represents the qualitat 
edge of dynamical systems theory to 



at 



lve 



Abelson's Bifurcation Interpreter [3]. 



uses MAPS to observe qualitative changes 
control parameters. These changes include: 



graphically displaying three-dimensional geo- 

of 3D polyhedral structures with different 

ics rendering tool, an engineer can visu- 

t, interactively edit the system to observe 

e useful portions of phase portraits with 

composite having the desired properties. 

the above scenario with the Phase Space 
esizes nonlinear controllers in phase spaces, 
able to synthesize a nonlinear control sys- 
sional and difficult to visualize, or on which 
e to obtain with traditional design tech- 
modeling of the phase-space stability regions 
our method. 



for 



od for automatically analyzing phase-space 

ems and have constructed MAPS to demon- 

a phase space, finds qualitatively different 

Ijrajectory flows with equivalence classes, and 
features. It employs deep domain knowl- 

r|ecognize the qualitative structures of phase 
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spaces. It computes a high-level descr|ipt 
bination of numerical, combinatorial, 
the phase-space structure with a relational 
of dynamical systems and the computational 
neering analysis enable us to build comput 
high-quality controllers for nonlinear 
MAPS to analyze a collection of ph 
control parameter values. 



lase 



ion of a dynamical system through a com- 

and geometric computations and represents 

graph. The qualitative representation 

formulation of the knowledge of engi- 

er programs that automatically design 

systems. In the next chapter we will use 

*. spaces each of which corresponds to fixed 
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-t^nase 



Phase 



Chapter 4 

Automatic 
Synthesis - 
Navigator 

4.1 Introduction 



This chapter develops Phase Space 
synthesis of nonlinear dynamical syst 
ically designs a controller for a nonlin 
control laws by synthesizing the desired 
and planning and navigating the syst 
space. 

The Phase Space Navigator relies 
control systems. It employs the MAPS 
tive phase-space structures characteri 
The control synthesis utilizes flow 
global control paths. We will present 
a paradigm for high-performance non 
trajectory planning and tracking will 
illustrated with an example of syntl 
column. 



pipes 



-Space Control 
Space 



Navigator, an autonomous system for control 
ems. The Phase Space Navigator automat- 
ear system in phase space. It generates 
phase-space flow "shapes" for the system 
em along good control trajectories in phase 

on the phase-space knowledge of dynamical 

5 program to extract and represent qualita- 

zing the qualitative aspects of the dynamics. 

to model phase spaces and to search for 

the novel idea of phase-space navigation as 

inear control design. Algorithms for control 

be described. The synthesis method will be 

g anti-buckling control laws for a steel 



lesizim 
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4.2 Automatic Control Synthesis in Phase Space 



The Phase Space Navigator syntr. 
of view in phase space. The control 
"steering" of the system trajectory, 
sired state by a control signal. 



lesizes 



a control system from a geometric point 

of a dynamical system is interpreted as the 

Emanating from some initial state, to the de- 



4.2.1 Overview of the Phase Space Navigator 



of 



The Phase Space Navigator consists 
tory generator, and a reference 
optimal paths from an initial state to 
sequence of path segments connected 
rameter changes. A brute-force, fine-£ 
would be prohibitively expensive. Hi 
trajectory flows provide a way to 
and search for global control paths 
information about the phase-space 
The trajectory follower tracks the 
deviations, and resynthesizes the 
changes significantly. 



4.2.2 Intelligent navigation in phase space 



The control objective for a stabi 
along which the physical system can b 
there afterwards under control. We 
physical plants to be controlled 
the plants cannot be ignored. When 
found for the desired state with 
paths have to be synthesized. G 
of the systems are far away from the 
advantage of the underlying dynami 
jectory both locally and globally and 



a global control path planner, a local trajec- 
trajecftory follower. The global path planner finds 
the goal state in phase space, consisting of a 
at intermediate points where the control pa- 
rain search in high-dimensional phase spaces 
gh-level descriptions of the phase space and 
eff ciently reason about phase-space structures 
The local trajectory generator uses the flow 
trajectories to produce smoothed trajectories, 
planned reference trajectory, reactively corrects 
reference trajectory if the dynamics of the system 



lizatlion problem is to synthesize a control path, 

e brought to the goal state and made to stay 

are particularly interested in cases in which 

operate in large regions where nonlinearities of 
no global stabilization control laws can be 

traditional control techniques, composite control 

eomjstrically, this is the case where initial states 

desired state. Phase Space Navigator takes 

:s of the phase-space flows to plan the tra- 

switches control at carefully planned time 
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instances and places in phase space 1 
planning requires relatively smaller control 
trol properties otherwise impossible 
It is a small, opportunistic dynamical 
of phase-space structures. As soon 
desired state, a local linear controller 



as 



4.2.3 Planning control paths with flow pipes 



pad 



blem 



gi 



The geometric modeling of a phase 
control planning and navigation feaslbl 
actions, the search for a control 
a reachability problem, i.e., the pro 
segments each of which is under a sin 
in Figure 4.1. This point-to-point pi 
pipe representation of phase spaces: t 
while, switch to a new control action, 
arrive at the goal. 

To make this approach comput 
dynamical system indexed by different 
set of trajectory flow pipes. These fl 
other and pasted together to form a 
graph is a directed graph where no 
are segments of flow pipes. The 
the goal state are nodes in the graph 
graph is weighed according to traveling 
a single value or a range of values, 
paths is formulated as a search for s 

Since a flow pipe models a bundL 
a small deformation of the reference 
This is useful as the synthesized 
accommodate uncertainties and noise 



low 



des 



graph 



With 



1 Variable-structure control [23] also corap 
our method differs fundamentally from the 
trajectory flows are modeled and utilized in 



This type of control via phase-space path 

authority and achieves the desired con- 

o obtain or difficult to manually synthesize. 

alteration based on the global knowledge 

the system enters the neighborhood of the 

stabilizes the system at the desired location. 



space with flow pipes makes the phase-space 
e. Given a discrete set of possible control 
ii from an initial state to a destination is 
of finding a sequence of connected path 
e control action, as schematically illustrated 
anning can be naturally executed in the flow- 
lie system can travel along one flow pipe for a 
jump onto another flow pipe, and eventually 

:ationally feasible, the phase portraits of the 

control actions are first parsed into a discrete 

pipes are then aggregated to intersect each 

graph, the flow-pipe graph. The flow-pipe 

are intersections of flow pipes and edges 

may contain cycles. The initial state and 

Each edge, a segment of a flow pipe, of the 

time, smoothness, etc. The weight can be 

this representation, the search for optimal 

hortest paths in the directed graph. 

of similar trajectories, it provides room for 

trajectory within the flow pipe at runtime. 

trajectory often has to be modified to 



reference 



oses phase spaces along switching surfaces. However, 
variable-structure control in how the phase-space 
the search for global control paths. 
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phase portrait 1 



phase portrait 2 



phase portrait 3 



Figure 4.1: Search for a control path fj-om an initial point to a goal point in a stack 
of phase portraits. 
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control 
objective 



Planning Module 

~l 



Optimality 
Constraints 



Nominal 
Model 
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Figure 4.4: Reference trajectory tracking 
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4.3.3 The autonomous control synthesis algorithms 

We consider a control system x' = f(x,u), where x G R n is the state variable 
and u £ R m is the control parameter. The initial state of the system is x s and 
the goal state for the system to reach is x g . The set of admissible values for the 
control parameter u is {u { }. The planning algorithm searches for paths from x s to 
x a , subject to certain design constraints, in the phase spaces indexed by control 
parameter values {uj. The tracking 
path. 

Algorithm 4.1 The trajectory plann 



algorithm follows the synthesized reference 



ling algorithm: 

1. Generating phase portraits and parsing them into flow pipes: 

Generate a collection of phase portraits indexed by {«,}. Parse each of the phase 
portraits into flow pipes 

2. Constructing flow-pipe graph: 

Aggregate flow pipes to form flow^pipe graph. Make x s and x g to be nodes in the 
graph. Weigh each edge of the graph with the cost of traversing the edge. 

3. Testing reachability of goal x t 

Cluster the flow-pipe graph into flow-pipe path components. If x s and x g are in the 
same flow-pipe path component, go to Step 6. Otherwise, the goal is unreachable 
with the given set of parameter values {«,-}; continue. 

Searching for partial paths from x s to x g : 

Compute the reachable set R Xs from x s and the reverse reachable set R x from x g in 
the flow-pipe graph. Denote the gap between R Xs and R Xg to be G = \\R Xt - R \\. 
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5. Tuning control parameter: 



Search for a value of the control 
sponding phase portrait has flow 



parameter u outside {ui] such that the corre- 
pipes that, when added to the flow-pipe graph, 
reduce the value of G. If G is zero, i.e., the gap between R Xs and R Xg is bridged, 
collect the flow-pipe paths from x s to x g and go to Step 7. Otherwise, repeat 
Step 5 until there are no more parameter values to search, in which case return 
the collection of partial flow-pipe paths from x s to x g . 



6. Finding flow-pipe paths from 

If looking for the shortest paths 



from x s to x g , run a standard shortest path 
algorithm on the graph. If searching for all paths p from x s to x g subject to the 
constraint C(p) < C, use the following algorithm 2 . 
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C( Pi ) + W(x it x)<(h; 
P'vi ~ {Pi * (xi,x)\x one edg^ reachable from Xi, x = x g , 

C(p t ) + W( Xi ,x)< 



par 



cy, 

tial— paths — \* partial— paths ~ Pi) U r pi ; 



(6b) 



-* global— P aths — * global— path 

(c) if Ppartial-paths / 0, goto 

(d) if Pghbal-paths = 0, goto Step 

(e) order paths in P g iobal- P aths 



X g l»CJ X q • 



paths — lA^'S} %s) J • 

ths, ending at X{, do 

: reachable from Xi, x =£ x g , 



P' . 

x Pi' 



; otherwise, continue. 
4; otherwise, continue, 
according to their costs and return. 



7. Generating trajectory segments: 

Select a representative trajectory segment from each flow pipe forming the flow- 
pipe path from x s to x g and paste them together at intersections of flow pipes 

8. Smoothing trajectories: 

Smooth each intersection of trajectory segments through trajectory homotopy de- 
formation. Order the smoothed trajectories with costs and return. 

Algorithm 4.2 The trajectory tracking algorithm: 



2 Note that in the discussion C(pt) is the (jost of a path p if W(x { , Xj ) is the weight of the edge 
(xi,Xj), and "*" is the path composition operator. 
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update the nominal model of the 



1. Sensing and estimating the physical system: 

Sense the state x of the physical system, estimate the model parameters, and 

system. If x is in the neighborhood of the goal 



x g , go to Step 5. If the nominal model changes substantially, signal and go to 
Step 4. Otherwise, continue. 

2. Computing control action: 

Based on the synthesized referente trajectory and the observation from Step 1, 
compute control tracking term u\ rack and correction term u correct . If u correct is 
unobtainable with a local linear correction, check alternative paths to the goal 
from the collection of suboptimal paths. If the goal is unreachable, go to Step 4. 
Otherwise, the control action is u\= u irack + u correct - continue. 

3. Generating control action: 

Tune the control parameter to u track and generate u correct with a local linear 
feedback controller. Drive actuators and go to Step 1. 

4. Resynthesizing reference trajectory: 

Call the planning module to resynthesize a reference trajectory and go to Step 2. 



5. Goal stabilizer: 

Stabilize the system at x g with a 
and noise. 



Figure 4.5 shows the flow chart of 
as input the system model, allowable 
and desired control objectives. It 
the form of a list of tuples: (time, 
The planning module has been 
the planning module synthesizes a 
The tracking algorithm described ab 
implementation is an immediate goal 



local linear controller, subject to uncertainties 



die planning algorithm. The algorithm takes 

control parameter values, design constraints, 

outputs a synthesized reference trajectory in 

switching point, control parameter value). 

implemented in Scheme. We will describe how 

law for stabilizing a buckling column. 

provides a conceptual framework, whose 

of future research. 



control 



ove 



4.3.4 Discussion of the synthesis algorithms 

On-line and off-line synthesis 

The planning module generates smoothed reference trajectories in the phase space. 



The constraints on response time and 



the availability of computational resources 
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Figure 4.5: The flow chart of the trajectory planning algorithm of the Phase Space 
Navigator. 
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dictate whether the computation of planning is done on-line or off-line. If the 
planning is done off-line, the synthesized plan is compiled into a table. At the run- 
time the controller performs a table lookup for a control action. On-line reference 
trajectory synthesis and tracking require substantial computational power. The 
Supercomputer Toolkit [4] provides a,n ideal platform for experimenting with the 
above control synthesis algorithm on real-time control applications. 

In Step 1 of Algorithm 4.1 for panning, the phase-space structure and flow 
descriptions can be either computed once for all the control parameter values or 
generated on demand in the search fo}- optimal paths. The latter case will be more 
suitable for on-line synthesis. 

Generalizations of the point-to-pioint planning 

We present the planning algorithm for a control design problem with one initial 
state and one goal state. The algorithm also applies to control systems of "one 
initial state/many goal states", "many initial states/one goal state", or "many 
initial states/many goal states". 

For a "one initial state/one goal sj^ate" or "one initial state/many goal states" 
problem, the Dijkstra's algorithm for the single source shortest paths in a directed 
graph runs in 0(V 2 ) for a graph with V vertices and E edges [11]. 

A "many initial states/one goal state" problem can be converted to the "one 
initial state/many goal states" problem by reversing the directions of all edges; and 
a "many initial states/many goal states" problem is solved by the Floyd- Warshall 
algorithm in 0(V 3 ) [11]. 

Discrete vs. continuous control parameter spaces 

A version of the synthesis algorithm is presented for the control parameter initially 
taking values from a finite discrete set. An example of such finite- valued control 
systems is the satellite attitude controller, which stabilizes the antenna direction 
of the satellite subject to disturbances by turning on or off high thrust jets. Other 
such examples are switching power regulators. 

The method also applies to control systems with continuous, multiple control 
parameter spaces. The program uniformly samples the continuous parameter space 
at discrete points and then applies the algorithm on these control points. A much 
better approach, however, would search for "land-mark points" that delimit dis- 
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tinct behaviors and then partition the parameter space into equivalent subspaces. 
For example, the Bifurcation Interpreter [3] can be employed to search for the 
bifurcation points defining qualitative changes in dynamics and decompose the 
parameter space into topologically equivalent subparts. 

Suboptimal control paths 

The algorithm finds optimal paths in the flow-pipe graph. However, the subop- 
timal paths can be useful when the optimal paths are no longer judged feasible. 
They can be stored in the table in addition to the optimal ones. A controller can 
opportunistically choose among available trajectory paths according to the desired 
properties at a control switching point. 

Comparison with dynamic programming 

To synthesize an optimal control path, dynamic programming discretizes a state 
space and conducts a fine-grain search in the discretization. The cost of the ex- 
haustive search could be prohibitive for large regions or in higher dimensions. In 
contrast, the Phase Space Navigator searches for control paths in a manageable 
set of flow pipes. It is possible to use dynamic programming within a flow pipe, 
once the global path has been established. 

4.4 An Example: Stabilizing a Buckling Column 

We illustrate the mechanism of the Phase Space Navigator with a control synthesis 
for stabilizing a buckling steel column. The buckling motion of the column has 
been extensively studied in nonlinear dynamics from a theoretical point of view 
and in structural engineering by practicing engineers. The columns are commonly 
used as strengthening elements in structures; for example, flexible space structures 
use columns in large operating regions. Study in nonlinear dynamics shows that a 
slender steel column can exhibit very complicated dynamical patterns under vari- 
ous operating conditions. Therefore, it is important to understand the dynamical 
behaviors of the columns and to devise ways of preventing them from failure. We 
will use the Phase Space Navigator to analyze the behaviors of the column and to 
synthesize a nonlinear controller for stabilizing the column. 
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Figure 4.6: Buckling of a thin elastic steel column due to axial end loads. 

4.4.1 The column model 

A thin elastic column, subject to axial end compressive loads, buckles around 
the principal axis, as shown in Figure 4.6. The nonlinearity is introduced by the 
nonlinear geometric stiffness of the column [31]. Stoker [47] gave a simplified model 
for the column subject to axial compressive force and viscous damping 

mx" + ex' -\- a^x + a 3 x 3 = 0, (4.1) 

where a x = A + C - 2P/1 and a 3 = B + D - P/l 3 . The state x is the characteristic 
measure of the column deflection from the principal axis and x' is the velocity. 
The column has mass m and length 21. The axial load is P, and the coefficient of 
viscous damping is c. The bending stiffness is modeled by a primary hard spring 
with restoring force Ax + Bx 3 and a secondary hard spring with restoring force 
Cx + Dx 3 . We rewrite equation (4.1) as a system of first-order equations 



(4.2) 



J X 1 — X-2 

\ X 2 = ^(-<*1^1 ~ a 3X 3 ~ CX 2 ), 

where x\ represents the deflection and x 2 represents the velocity. 

The system (4.2) describes the buckling motion of the column and represents 
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Figure 4.7: The buckling column: the phase space of a buckling column showing 
the stability boundaries and connecting trajectories. The horizontal axis x is the 
characteristic measure of the column displacement from its principal axis and the 
vertical axis x' is the velocity. 

only a single mode of vibration. For a long and slender column, vibrations are 
observed to occur primarily in the first mode [32]. 

When the axial load P is less than the critical load P cr iticah the column oscillates 
around the principal axis and returns to the vertical state. It has only one stable 
state corresponding to the attractor at the origin of the phase space. Under a 
heavier load, i.e., P > Pcriticah the column buckles to one side or the other. The 
phase space has a saddle at the origin and two attractors symmetrically arranged 
about the saddle; see Figure 4.7. The positions of the two buckled states depend 
on the external load P, the stiffness of the column, and the length /. The larger 
the load is, the farther away the buckled states are from the principal axis. The 
column breaks when the buckling exceeds certain limit. As P surpasses P cr iticah 
the column undergoes a pitchfork bifurcation of equilibria [18]: the attractor at 
the origin gives birth to a saddle at the origin and two attractors on two sides. 
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4.4.2 Extracting and representing qualitative phase-space 
structure of the buckling column 

MAPS automatically analyzes the column model in the phase space and extracts 
and represents the qualitative phase-space structure [54]. For parameter values 
ai/m = —2.0, a 3 /m = 1.0, and c/m = 0.2 and phase-space region —3.0 < x\ < 3.0 
and —4.0 < x 2 < 4.0, the program reports the following findings and represents 
them internally in a relational graph: 

<equilibrium-points : 
equilibrium 1. (attractor at (1.41 0.)) 
equilibrium 2. (saddle at (0. 0.)) 
equilibrium 3. (attractor at (-1.41 0.))> 

<trajectories: 
<boundary-traj ectories : 

trajectory 1. (from *infinity* to (0. 0.)) 

trajectory 2. (from *infinity* to (0. 0.))> 
<connecting-traj ectories: 

trajectory 3. (from (0. 0.) to (-1.41 0.)) 

trajectory 4. (from (0. 0.) to (1.41 0.))» 

<stability-regions : 
stability-region 1. 

attractor at *infinity* 

stability-boundary: () 

connecting-trajectories: () 
stability-region 2. 

attractor at (1.41 0.) 

stability-boundary: (trajectory 2. trajectory 1.) 

connecting-trajectories: (trajectory 4) 
stability-region 3. 

attractor at (-1.41 0.) 

stability-boundary: (trajectory 2. trajectory 1.) 

connecting-trajectories: (trajectory 3)> 

The program finds two attractors at (1.41,0.0) and (-1.41,0.0) and a saddle at 
the origin. It generates a high-level description of the phase-space geometry: two 
banded stability regions associated with the two attractors, separated by the sta- 
ble trajectories of the saddle at the origin. Figure 4.7 shows stability boundaries 



84 




Figure 4.8: The buckling column: the flow pipe leading to the attractor on the 
left. 

and connecting trajectories of the two stability regions. Based on the geomet- 
ric phase-space representation, the phase space is further modeled with two flow 
pipes formed by aggregating geometric pieces, using the flow-pipe construction 
algorithm. The pipe boundaries consist of the separatrices of the two stability 
regions that approach the saddle and of the trajectories that connect equilibria, as 
described by the program in the following: 

<f low-pipes : 
flow-pipe 1. (from *infinity* to (-1.41 0.)) 

boundary: (trajectory 2. trajectory 1. trajectory 3.) 
flow-pipe 2. (from *infinity* to (1.41 0.)) 

boundary: (trajectory 2. trajectory 1. trajectory 4.)> 

Figure 4.8 shows the flow pipe that ends at the left-hand attractor. 



4.4.3 Synthesizing control laws for stabilizing the column 

We want to stabilize the column at its vertical state to support the axial end 
loads and to prevent it from breaking. Under sufficiently heavy load, the buckling 
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motion of the column leads the column to one of the buckled states and, when the 
states are far away from the principal axis, induces the failure of the column. 

We shall focus on global navigation in phase space that synthesizes global 
reference trajectories leading to the desired goal state. Since the design of local 
linear controllers is relatively well understood, we shall not discuss the trade-offs 
among different designs for linear controllers and shall choose, for the purpose of 
demonstration, a simple linear feedback design. Local controllable regions of such 
linear controllers are quantified, given available control strength, and are used in 
constraining the design of global control paths. 

The goal state of the control is the unbuckled state — the saddle at the origin 
of the phase space that does not have a stability region. We want to synthesize 
a non-zero stability region for the goal and maximize the region. The controlled 
column is of the form 

{x[ = x 2 
X 2 = - l (- a l X l ~ Wi - CX 2 ) + IX, 

where u is the control. In the model, mu has the same dimension as the force. 

The Phase Space Navigator automatically synthesizes a global trajectory from 
an initial state x s to the goal state x g . We consider two cases for the initial state 

x s . 

Control design I: stabilizing the buckling motion 

The column is initially buckling with sufficient velocity. The initial state x s is far 
away from the saddle x g in phase space in this case. The control parameter u takes 
its values from a small range around zero. 

The Phase Space Navigator constructs a flow-pipe graph from phase spaces of 
u = and u =*local-control*, as shown in Figure 4.9. The graph shows the 
case when the initial state x s is in the flow pipe ending at the left attractor as 
shown in Figure 4.8. Denote by U ctl the local controllable region at the goal x g . 
The two flow pipes of u = intersect U ct i at m and u 2 respectively. The edges u{x: 



3 We use this example for the purpose of illustrating the Phase Space Navigator. The sys- 
tem (4.3) is actually feedback linearizable. A feedback linearization would cancel the nonlinearity 
of the vector field in x 2 direction. In contrast, our design restrains the control to be less than 
10% of the vector-field strength in the first case and less than half of the vector-field strength in 
the second case. 
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left-attractor (-1.41, 0) 




right-attractor (1.41, 0) 



infinity 5 



Figure 4.9: The flow-pipe graph for the buckling column. 

and u^T g represent flows within U ct i produced by the local controller. The paths 
(^infinity* -> x s -> m -> left-attractor) and (^infinity* -► u 2 -» right-attractor) 
represent the two flow pipes of the phase space of u = 0. 

The Phase Space Navigator finds a simple path from x s to x g in the graph, 
consisting of edges x^T[ and u{x~ g . Then the program synthesizes an individual 
trajectory connecting x s and x g from the flow-pipe path (xJT[, W^)- To construct 
this desired trajectory, the program deforms an uncontrolled trajectory, emanat- 
ing from x 8 , of the flow pipe represented by the edge x^u~[ so that the deformed 
trajectory enters U ctl . The region U ct i is first projected backwards through the flow 
pipes to form a goal projection. Figure 4.10 shows the goal projection from U ctl — 
two thin pipes illustrated in thick solid lines. Then the uncontrolled trajectory is 
deformed towards the goal projection. The programmed deformation is designed 
to push the trajectory towards the nearest goal projection in the controllable di- 
rection x 2 . Since the control is more effective when the direction of the vector 
field is relatively orthogonal to the controllable direction x 2 , switching points are 
inserted to turn off the controller when the angle between the two directions are 
below some threshold. 

Consider for example the case when the column is initially at the state (-1, -3). 
The control design is specified as: 
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Figure 4.10: The goal projection and the deformation of the trajectory. 



control_type : 
goal_state: 
initial_state: 
range_of .control : 



point_to_point 
(0.0, 0.0) 
(-1.0, -3.0) 
u € [-0.2,0.2] 



The Phase Space Navigator designed a control law that brings the column back 
to the unbent state. The control law is represented as a list of tuples, each of which 
specifies the time, state, and control value for each switching of control: 

((time 0.) (switching-state #(-1 -3)) (control .2)) 

((time .284) (switching-state #(-1.82 -2.71)) (control 0.)) 

((time 1.06) (switching-state #(-1.86 2.49)) (control -.2)) 

((time 2.71) (switching-state #(1.35 1.82)) (control 0.)) 

((time 6.76) (switching-state #(-.0023 -.0692)) (control *local-control*) ) 

Figure 4. 11 (a) shows the synthesized reference trajectory originating at (-1,-3), 
The circles indicate places where the control of deformation switches. Each seg- 
ment of the reference trajectory delimited by the switching points is under interval- 
constant control, as specified by (71, £72, (73, (74, or a local linear control law as 
the trajectory is in the vicinity of the goal. The global portion of the reference 
trajectory is pushed towards the goal projection with small deformation that is less 



than 10% of the vector field strength at any state or 0.2 in the normalized unit, 
whichever is smaller. The position x and velocity v = x' of the controlled column 
are plotted against the time t in Figure 4.11(b) and Figure 4.12(c), respectively. 
The control u is shown in Figure 4.12(d). 

When the goal is not reachable with the current deformation, the synthesis 
algorithm tries again with increased control strength if possible. The synthesized 
trajectories could be further optimized with variational techniques on the collection 
of trajectories within the flow-pipe segments [6]. 

Control design II: restoring from the buckled state 

The column is initially near the buckled state (-1.41,0.0). The control objective 
here is to pull the column out of the buckled state and to bring it close to the 
unbuckled state. Since the flow pipe containing the initial state does not intersect 
U c ti in the down stream, a different control strategy must be employed. The 
program uses the following strategy to synthesize a control path. 

Assume the initial state of the column is (-1.5,0.0). The control objective is: 

control_type: point_to_point 

goal.state: (0.0, 0.0) 

initial_state: (-1.5, 0.0) 

range_of .control : ue [-1.0, 1.0] 

The uncontrolled trajectory from the initial state would spirally approach the 
attractor in the clockwise fashion. The control is exerted in such a way as to swing 
the trajectory away from the buckled state to approach the goal projection of the 
saddle. When the trajectory intersects the goal projection, the control is switched 
off so that the system slides along the uncontrolled trajectory. As soon as the 
system enters U ct i, the linear controller is switched in. The control strength is less 
than half of the vector field strength or 1.0 in the normalized unit, whichever is 
smaller. The synthesized control law is specified as: 

((time 0.) (switching-state #(-1.5 0)) (control -.187)) 
((time .001) (switching-state #(-1.5 .000187)) (control .187)) 
((time 1.65) (switching-state #(-1.24 -.0011)) (control -.289)) 
((time 3.19) (switching-state #(-1.66 .0153)) (control .638)) 
((time 5.92) (switching-state #(-.527 -.00186)) (control -.454)) 
((time 7.74) (switching-state #(-2.02 .0473)) (control 1.)) 
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(a) 




(b) 



Figure 4.11: The synthesized control law that stabilizes the buckling column: (a) 
the reference trajectory that leads to the unbuckled state corresponding to the 
saddle at the origin. The column is initially buckling with sufficient velocity; (b) 
the position x of the column plotted against time t. (to be continued) 
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(c) 




(d) 



Figure 4.12: (cont'd) The synthesized control law that stabilizes the buckling col- 
umn: the velocity v(= x') of the column and the control signal u for stabilizing 
the column are plotted against time t in (c) and (d), respectively. 
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((time 8.08) (switching-state #(-1.75 1.47)) (control 0.)) 

((time 11.1) (switching-state #(-.049 .00163)) (control *local-control*) ) 

The corresponding reference trajectory is shown in Figure 4.13(a). The control 
parameter changes at places marked by circles. The reference trajectory consists of 
eight trajectory segments labeled by interval-constant control (71, (72, ..., (77, and 
a local control law. The first segment (£71) and the last segment (local control) are 
very short in length and thus are invisible in the figure. The position x and velocity 
v = x' of the controlled column are plotted against the time t in Figure 4.13(b) 
and Figure 4.14(c), respectively. The control u is shown in Figure 4.14(d). 

Control design III: local control near the goal 

A simple local linear controller of the form u local _ Hnear = -k lXl - k 2 x 2 is used to 
stabilize the system around the goal (see [55] for details). Given the maximum 
control strength 0.2 for the linear controller, we choose an over-conservative region 
for U ct i determined by 

U c ti = {(x u x 2 )\ 10.835a;! + 0.550a; 2 1 < 0.04, |0.797ar! — 0.605x 2 | < 0.04}. (4.4) 
This U ct i is used in both of the above cases. 

4.4.4 The phase-space modeling makes the global naviga- 
tion possible 

The qualitative description of the phase-space structure and the geometric mod- 
eling of the trajectory flows provide a "map" for navigating system trajectory to 
the goal in phase space. 

The directed graph constructed from the flow pipes is used in searching for 
global paths and in determining whether the goal is reachable. The flow pipes 
also make it possible to characterize the more microscopic spatial and temporal 
relations between the current state and the goal state. For example, the collection 
of simplices near the goal forms a neighborhood of the goal that determines whether 
or not the trajectory has missed the goal in the flow pipe. The deformation of 
global path segments is constrained by reasoning about the spatial relation with 
flow pipes. 
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(a) 




(b) 



Figure 4.13: The synthesized control law for restoring the column from the buckled 
state: (a) the reference trajectory that swings the column out of the buckled state; 
(b) the position x of the column plotted against time t. (to be continued) 



93 




(c) 




(d) 



Figure 4.14: (cont'd) The synthesized control law for restoring the column from 
the buckled state: the velocity v(= x') of the column and the control signal u for 
controlling the column are plotted against time t in (b) and (c), respectively. 
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4.4.5 Other control problems 

The stabilization for a buckling column is closely related to the pole-balancing 
problem, where the control objective is to balance a pole at the vertically upward 
state. In the phase space of the pole model, the saddle is the upward state of 
the pole and the two attractors on two sides of the saddle correspond to the same 
downward state of the pole. 

The control objective for pole-balancing is to bring the pole to the upward 
state from any other positions, even when the pole is initially hanging downward. 
The problem has been studied as a standard testbed for many control designs, 
most of which focus primarily on linear feedback control in linear regimes. With 
our phase-space navigation strategy, a globally stable reference trajectory can be 
synthesized, similar to the one for the buckling column discussed earlier. 

The optimal design of such trajectories has many practical implications. For 
example, the efficiency of cargo handling work at shipyards depends largely on the 
operation of overhead cranes [44]. Figure 4.15 illustrates such a crane driven by a 
trolley drive motor and a hoist motor. The bottleneck of the ship unloading oper- 
ation is the transfer of cargo from a ship to waiting trucks; therefore, minimizing 
this transfer time brings about a large cost saving. The planar motion of the load 
hanging at the moving trolley can be modeled as a swinging pendulum with a mov- 
ing point of attachment and a hoisting cable. The trajectories of such overhead 
cranes could be optimally designed with the Phase Space Navigator. This would 
be an important step towards the full automation of cargo handling work without 
a crane operator. 

4.5 Summary of the Chapter 

We have developed an autonomous control synthesis method, the Phase Space 
Navigator, for synthesizing controllers for nonlinear systems in the phase spaces. 
We have discussed primarily finite- valued control systems as examples to illustrate 
the method and have noted that the method also applies to control systems with 
continuous parameter spaces. 

The Phase Space Navigator synthesizes global reference trajectories using knowl- 
edge of phase-space structures provided by the modeling and analysis program. A 
phase space is modeled with flow pipes that are collections of trajectories hav- 
ing the same qualitative behaviors. The flow pipes provide a way to efficiently 
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Figure 4.15: The mechanical model for an overhead crane unloading cargo from 
ships, reproduced from [Sakawa&Shindo, 1982]. The crane is equipped with a 
trolley drive motor and a hoist motor. The planar motion of the load is modeled 
as a swinging pendulum hanging at the moving trolley. 

search for and reason about global structures of phase spaces. The global nonlin- 
ear control synthesis becomes a graph search problem with this representation of 
the phase space. We have shown how global reference trajectories can be automat- 
ically designed with the example of stabilizing a buckling column. With the novel 
idea of grouping an infinite number of trajectories into a manageable collection of 
flow pipes, the difficult task of synthesizing a nonlinear controller is formulated 
as a computational problem that requires a combination of substantial numerical, 
symbolic, and combinatorial computations and spatial reasoning techniques. 

The leverage of the automatic synthesis method comes from applications whose 
operating regions are grossly nonlinear and on which high-performance global con- 
trollers are impossible to synthesize with traditional techniques. In the next chap- 
ter, we will demonstrate the application of the Phase Space Navigator to the design 
of a maglev controller. Other potential engineering applications include large flex- 
ible space structures, robot manipulator planning, satellite attitude control, and 
switching power regulators. 
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Chapter 5 

An Application: Design of a 
Maglev Controller 



This chapter studies a magnetic levitation system and applies the machineries of 
the Control Engineer's Workbench developed in Chapters 3 and 4 to the design of 
a stabilizing controller for the system. 

The Workbench has helped synthesize a global, nonlinear controller for the 
nominal model of a German maglev system 1 . We describe the systematic state- 
space design method for determining the global switching points of the controller. 
The synthesized control system can stabilize the maglev vehicle with large initial 
displacements from an equilibrium. The simulation shows that our nonlinear con- 
troller possesses an operating region more than twenty times larger than that of 
the classical linear feedback design for the same system. 

5.1 Introduction 

Magnetically levitated trains provide a high speed, very low friction alternative to 
conventional trains with steel wheels on steel rails. Several experimental maglev 
systems in Germany and Japan have demonstrated that this mode of transporta- 
tion can profitably compete with air travel. More importantly, maglev transporta- 
tion can ease traffic congestion and save energy [14, 28, 50]. 



! The maglev stabilization problem was suggested to the author by Professor Richard Thornton 
of Electrical Engineering at MIT, who is leading the MIT Maglev Consortium. The result of the 
work described here appears in a paper [56] to be presented at the 31st IEEE Conference on 
Decision and Control, December 1992. 
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Maglev transportation uses magnetic levitation and electromagnetic propul- 
sion to provide contactless vehicle movement. There are two basic types of mag- 
netic levitation: electromagnetic suspension (EMS) and electrodynamic suspension 
(EDS). In EMS, the guideway attracts the electromagnets of the vehicle that wraps 
around the guideway. The attracting force suspends the vehicle about one cen- 
timeter above the guideway. In contrast, the EDS system uses repulsive force, 
induced by the magnets on the vehicle, to lift the vehicle. 

5.2 The Maglev Model 

An attractive system such as the EMS system is inherently unstable. We consider 
the control design for stabilizing an EMS-mode train traveling on a guideway— a 
simplified model for the German Transrapid experimental system. The Transrapid 
system is schematically shown in Figure 5.1. It uses attractive magnetic forces to 
counterbalance gravitational forces. 

The state equations for the magnetically levitated vehicle and the guideway are 
described by 

dx z(Vj-Rx) , xy 

dt L z ' z 

f=^-^ (5.1) 

where the state variables x, y, and z represent coil current in the magnet, vertical 
velocity of the vehicle, and vertical gap between the guideway and the vehicle, 
respectively. The control parameter is the coil input voltage V l . The other param- 
eters are the mass of the vehicle m, the coil resistance R, the coil inductance L 
and the vertical gap z at the equilibrium, and the gravitational acceleration g. 
Details of the derivation of the model are discussed in [49]. The nonlinearities of 
the system come from the nonlinear inductance due to the geometry of the magnet 
and the inverse square magnetic force law. 

The system has one equilibrium state at which the magnetic force exactly 
counterbalances the force due to gravity and the vehicle has no vertical velocity and 
acceleration. However, the equilibrium is a saddle node which is not stable. The 
control objective, therefore, is to stabilize the vehicle traveling down the guideway 
and maintain a constant distance between the vehicle and the guideway despite 
any roughness in the guideway. The available control input is the coil input voltage 
Vi in the model (5.1). We further assume that V t is produced by a buck converter 
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Figure 5.1: EMS maglev system for high-speed ground transportation, represent- 
ing a simplified drawing of the German Transrapid design, (a) Electromagnetic 
suspension (from [MTAC Report 1989]); (b) Detail of a suspension magnet, super- 
imposed on the field distribution (from [Eastham 1989]). 
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capable of delivering any voltage from to 300 volts. 

A linear control design for the maglev system described in [49] uses the pole- 
placement method. The system is first linearized around the equilibrium. The 
linearized system has unstable poles, i.e., the poles in the right-half of s-plane. 
A linear feedback is introduced to move the poles to the desired locations in the 
left-half of the s-plane. Such a control design can bring the system back to the 
equilibrium with an initial displacement of up to 0.2mm from the equilibrium. The 
linear controller saturates at the beginning for larger initial displacements. This 
is because the linearized model no longer approximates the original system well 
in regions far away from the equilibrium. A global, nonlinear control law such as 
a bang-bang control that respects the nonlinearity of the system must therefore 
precede the linear feedback control. However, the real challenge for the nonlinear 
design is to determine the global control law specifying, for instance, the switching 
points. 

5.3 State-Space Control Trajectory Design 

We will describe a nonlinear control design — a switching-mode control — in state 
space for the maglev system with large initial displacements from the equilib- 
rium. We will show that this controller can be automatically designed with the 
Workbench comprising MAPS and Phase Space Navigator developed earlier. The 
nonlinear controller brings the system to the vicinity of the equilibrium and then 
switches to the linear controller. 

For the purpose of demonstration, we assume that the vehicle is displaced 
from the equilibrium in the direction further away from the guideway. We will 
concentrate on the global design of the control reference trajectories and assume 
that a linear feedback design is available as soon as the system enters the capture 
region of the linear controller. 

5.3.1 Modeling state-space geometry 

The global control law is designed by analyzing and modeling the state-space 
geometry of the system. The Workbench explores the state space of the system 
and characterizes the state space with stability regions and trajectory flow pipes. 
It composes the state spaces for diiferent control parameter values and uses flow 
pipes to synthesize a composite state space. 

100 



The state variables x,y, and z in the model are scaled by 1, 10 3 , and 2 X 10 4 , 
respectively. The parameters of model are assumed to be: L = O.lh, z = 0.01m, 
R = lH, m = 10000%, and g = 9.8m/ sec 2 , typical of a large vehicle lift magnet. 
Assume the power supply delivers 140 volts, i.e., V { = 140, at the equilibrium. The 
Workbench explores the state space in a region bounded by the box {(x,y,z)\x G 
[0,400], y € [-300, 350], ^ G [0,600]} and finds the following equilibrium point: 

saddle: #(140. 0. 200.) 

eigenvalues: -17. 004+22. 963i 

-17. 004-22. 963i 

24.007 
eigenvectors: #(.23604 .97174 0) 

#(.51331 -.55588 .65384) 

#(.30157 .73255 .61027) 

With the information about the stable eigenvectors of the saddle, the Work- 
bench computes the stable manifold of the saddle, a two-dimensional surface. The 
Workbench generates a set of trajectories evenly populating the stable manifold to 
approximate the surface. The trajectories are obtained by backward integrations 
from initial points in a small neighborhood of the saddle. This neighborhood lies 
within the plane spanned by the stable eigenvectors of the saddle. 

Figure 5.2 shows the trajectories on the stable and unstable manifolds of 
the saddle in the ^-projection of the state space. The stable manifold is two- 
dimensional and the unstable one is one-dimensional. The stable manifold sep- 
arates the state space into two halves: trajectories in the upper-half approach 
z -> oo along one of the unstable trajectories, corresponding to the case in which 
the vehicle falls off the rail; and trajectories in the lower-half approach z = plane 
along the other unstable trajectory, corresponding to the case in which the train 
collides with the rail. 

5.3.2 Synthesizing a global stabilization law 

For an initial displacement above or below the equilibrium, the uncontrolled system 
will follow either a trajectory traveling upwards with increasing z and leaving the 
bounding box or a trajectory traveling downwards and hitting the z = plane. To 
stabilize the system at the equilibrium, it is necessary to synthesize a new vector 
field on both sides of the stable manifold so that trajectories travel towards the 
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Figure 5.2: Stable and unstable manifolds of the saddle for V t = 140 (yz- 
projection). 

stable manifold of the saddle in the new vector field. We consider only the top-half 
here. 

The Workbench first explores the state space of the model for different values 
of Vi and concludes that the larger the V- is, the further away the stable manifold 
is from the z = plane. For V t = v > 140, the region sandwiched by the stable 
manifold of Vi = 140 and that of V = v has the desired property— the vector field 
of Vi = v in this region is pointed towards the stable manifold of V = 140. When 
v = 300, the region is maximized. 

Similarly, the Workbench finds that the model with K = 300 has a saddle node 
at (300., 0., 428.57). The stable manifold of the saddle has a similar structure as 
that of Vi \ = 140 case, but with larger z coordinate. The ^-projection of the stable 
manifold and unstable trajectories for Vi = 300 is shown in Figure 5.3. Appendix A 
contains the information about the saddle for K = 140 and Vi = 300. 

Figure 5.4 shows the t/z-projection of the sandwiched region discussed above. 
The region is bounded by three pieces of triangulated surfaces, the top boundary 
shown in Figure 5.5(a), the side one in Figure 5.5(b), and the bottom one in 
Figure 5.5(c), respectively. The top boundary represents the stable manifold of 
the saddle for Vi = 300 and the bottom one for K = 140. The trajectories of 
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Figure 5.3: Stable and unstable manifolds of the saddle for V { = 300 (yz- 
projection). 

Vi = 300 flow into the region from the side boundary in Figure 5.5(b) and leave 
the region at the bottom boundary in Figure 5.5(c). There is no flow across the 
top boundary in Figure 5.5(a). 

The Workbench determines that the above region consists of two trajectory 
flows: the one for V- = 300 that flows into the region on the side boundary and 
pierces through the bottom boundary, and the one for V t = 140 on the bottom 
boundary that approaches the desired equilibrium in the limit. With the flow-pipe 
characterization of the state-space trajectory flows, the Workbench searches for 
control trajectories in this set of flow pipes and finds a sequence of flow pipes that 
lead to the desired goal: the composite of the trajectory flow for V- = 300 and the 
flow for Vi = 140, glued together at the stable manifold of V- = 140 represented by 
the bottom boundary of Figure 5.5(c). As a result, all the trajectories of V { = 300 
within the region can be brought to the equilibrium by switching to V- = 140 
as soon as the trajectories hit the bottom boundary. We call this region the 
controllable region for the system, as shown in Figure 5.4. 

Since the bottom boundary is an approximation to the true stable manifold, 
the trajectories of V- = 140 on the boundary can only get close to the desired 
equilibrium. The closeness depends on the quality of the manifold approximation. 
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Figure 5.4: The sandwiched region in ^-projection. 
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Figure 5.5: The boundaries of the sandwiched region in T/z-projection: (a) top 
boundary; (b) side boundary; (c) bottom boundary. 
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As the trajectories enter a small neighborhood of the equilibrium, we use a linear 
feedback controller, such as the one discussed in [49], to stabilize the system at 
the equilibrium. Figure 5.6 shows the synthesized control reference trajectories 
originating at different initial displacements from the equilibrium: 1mm, 4mm, 
4.5mm, and 5mm. The controller is able to recover from the first three initial 
points that are within the region. The last point is outside the region and thus 
uncontrollable; the current in the magnet can not build up fast enough to keep 
up with the ever-increasing airgap. The complete print-out of the synthesized 
reference trajectories showing the switching points is in Appendix B. 

For example, the control law for the initial displacement of 1mm is specified 



as: 



((time 0.) (switching-state #(140. 0. 220.)) (control 300.)) 

((time .0134) (switching-state #(163. 4.44 221.)) (control 140.)) 

((time .127) (switching-state #(139. .0789 202.)) (control *local-control*)) 

This control law specifies the time instance for each switching, the switching 
state, and the corresponding control value during the following time interval. The 
corresponding reference trajectory consists of trajectory segments delimited by 
these switching states. 

5.3.3 Evaluating the control design 

The synthesized global control law is a switching-mode one that changes the control 
parameter at the switching surface— the stable manifold of V { = 140. It is able to 
bring trajectories originating from any states within the controllable region to a 
local neighborhood of the saddle. 

The responses of the controller with respect to the four different initial displace- 
ments are shown in Figures 5.7 and 5.8. The vertical axis of each graph represents 
state variables x, y, and z as in the maglev model (5.1), one for each curve, and the 
horizontal one represents the time. For all the controllable initial displacements, 
the controller is able to bring the system back to the equilibrium with errors less 
than 0.2mm in displacement— a distance within the capture range of the linear 
feedback controller. 

By exploring the state-space geometries of the maglev system, the Workbench 
is able to automatically determine the switching points for the global controller. 
The linear feedback controller can recover from only displacements of less than 
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Figure 5.6: The synthesized control reference trajectories originating from four 
different initial states, together with the controllable region: (a) ^-projection; (b) 
zrc-projection. 
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0.2mm. The global controller has significantly enlarged the operating region of the 
linear controller. With the geometric representation of the controllable region in 
state space, the Workbench precisely determined that the maximum recoverable 
displacement is 4.55mm. Our simulation geometrically explained the observation 
in [49] that the vehicle would fall off the rails with a displacement of 5mm or larger. 
Many issues remain to be addressed in order to make the control design prac- 
tical. Since our control law is designed with the nominal model of the maglev 
system, the effect of uncertainties in the model and of noise in the environment on 
the design needs to be studied in future research. The design can also be optimized 
with respect to response time. 

5.3.4 Visualizing the design 

The synthesized controllable region is modeled with a polyhedral structure. This 
structure can be presented to engineers in a visual way. J. Choi has implemented 
a graphic rendering program for this purpose. Figure 5.9 shows a picture of the 
graphically rendered controllable region. The graphical presentation facilitates in- 
teractive, incremental modifications to the design in terms of the geometry. For 
example, one might want to directly manipulate the controllable region by tuning 
knobs corresponding to parameters. The geometric representation helps the de- 
signers to develop intuitions and to visually explore the effects of certain design 
choices. 

5.3.5 Implementation of the controller 

The control design described above has been computationally simulated only. How 
will this design be implemented on the real system and used in real time? 

The control law specifying the switching surfaces in state space can be com- 
piled into a table. The control execution will be a table lookup and a geometric 
inequality test. At each step of the execution, the state of the system is sensed 
and checked against the switching surface. If the state is on the switching surface, 
the corresponding control value for the next time interval is read from the table 
and applied to the physical system. The implementation does not have to be very 
different from that for a dynamic programming one. 
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Figure 5.7: Simulation of the nonlinear control design for different initial displ 
ments. (to be continued) 
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(d) dz = 5mm 

Figure 5.8: (cont'd) Simulation of the nonlinear control design for different initial 
displacements. 
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Figure 5.9: The controllable region for the maglev control design, rendered with 
the light source straight on. Hidden faces and lines are removed. 
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5.4 Summary of the Chapter 

We have described the automatic design of a high-quality controller for stabilizing 
an EMS-mode maglev vehicle traveling above a guideway with the Workbench. 
The synthesized control reference trajectory consists of a sequence of trajectory 
segments, connected at intermediate points where the control voltage switches. At 
run-time, the controller will cause the system to track the reference trajectory and 
reactively correct local deviation from the desired trajectory. We have illustrated 
the state-space method for designing the global switching points of the nonlinear 
controller. The simulation showed that our design allows the maglev train to oper- 
ate with much larger disturbances in the airgap than the classical linear feedback 
design does. 
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Chapter 6 



Related Work 



This chapter surveys relevant work in areas of simulation and analysis for complex 
dynamical systems and nonlinear control synthesis methods. Whenever possible, 
the Workbench will be compared with related systems elsewhere. 

6.1 Qualitative Analysis of Dynamics 

Qualitative analysis of dynamics is concerned with deriving symbolic descriptions 
and constraints for the dynamics of a physical system from simulations. The MAPS 
program of Chapter 3 falls into this category. Abelson and Sussman described a 
collection of computer programs that automatically analyze dynamical systems at 
the level of expert dynamicists [1, 2]. 

Yip has constructed a program, KAM, for automatically analyzing Hamilto- 
nian systems with two degrees of freedom in two-dimensional phase sections [52]. 
The program uses techniques from computer vision to cluster trajectory points in 
phase sections and classifies phase portraits into meaningful categories. MAPS 
works in the domain of dissipative dynamical systems; it analyzes the systems in 
their continuous phase spaces regardless of their dimensions. In contrast, KAM 
applies to Hamiltonian maps on two-dimensional phase sections. MAPS constructs 
a simplicial representation for the phase-space qualitative features based on geo- 
metric pieces, as required by the task of control design. KAM uses a point-set 
representation instead. 

Sacks' Poincare program analyzes planar systems through a partition algo- 
rithm on phase spaces and a bifurcation analysis on one parameter [43]. The 
partition algorithm determines flow patterns in the plane using the properties of 
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two-dimensional flows. These properties are based on the Jordan Curve Theorem 
and the fact that trajectories are one-dimensional curves that do not intersect. 
However, the special properties about two-dimensional flows do not generalize to 
higher dimensions. Recently, Sacks has been exploring the computation of two- 
dimensional manifolds in three dimensions using a set of trajectories, similar to 
that of MAPS. MAPS does not analyze bifurcations. MAPS differs from Sacks' 
partition algorithm in that MAPS is able to analyze phase spaces of any dimen- 
sion, based on a general theoretical result on dynamical systems. In addition to its 
reasoning about the phase-space flow patterns, MAPS builds a workable represen- 
tation in the form of a relational graph characterizing the spatial arrangement of 
phase-space objects. The representation contains information about the geometries 
of stability regions. 

Eisenberg's Kineticist's Workbench is a program for simulating and explain- 
ing chemical reaction mechanisms [15]. It analyzes a chemical mechanism with 
numerical integrations and graph algorithms and interprets the simulation results 
in an episode data structure. The program does not use phase-space properties 
of chemical dynamical systems. Nishida et al. described a program for analyz- 
ing second-order nonlinear systems in phase space [34]. The program uses a flow 
grammar to specify possible patterns of solution curves in the phase plane and is 
restricted to two-dimensional flows. 

Hsu developed a discrete, cell-based method called cell-to-cell mapping for ap- 
proximating state spaces [21]. A continuous state space is discretized into regular 
cells forming a cell space. The associated map of a system becomes a cell-to-cell 
map which maps one cell to another cell. The cell-to-cell mapping method ap- 
proximates the stability region of an attracting cell with a collection of cells that 
eventually map to that cell. MAPS does not take this uniform approach. Instead, 
it partitions a phase space into discrete regions according to trajectory flow pat- 
terns. The partition respects the geometry of the phase space. Because of the 
hierarchical (top-down) nature of Hsu's method, it can be integrated with MAPS 
in refining an approximation for a stability region. 

Much of the difference between MAPS and the above cited work is due to the 
fact that the analysis MAPS performs is meant for control design. The representa- 
tion in MAPS makes explicit the effects of certain design choices and is meaningful 
to control designers. For example, the geometries of stability regions are modeled 
m this representation; the sizes of the regions are directly linked to certain design 
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parameters. 

6.2 Engineering Stability Analysis 

Another popular class of methods for engineering stability analysis is the so-called 
Liapunov-function based constructions for stability regions (domains of attraction). 
Our stability analysis of Chapter 3 does not use Liapunov functions; it is based on 
the stability-boundary characterization by Chang et al. 

The class of Liapunov-function based approaches arises from the Liapunov the- 
ory of stability. It estimates domains of attraction by constructing suitable Lia- 
punov functions. A Liapunov function is a family of closed surfaces such that a 
system trajectory remains in the region bounded by a surface after the trajectory 
enters. The Liapunov function is difficult to construct. Numerous algorithms have 
been proposed for this construction. Reference [17] gives a survey of these algo- 
rithms. Most of the algorithms assume special forms of the Liapunov functions 
and determine the unspecified coefficients in the functions. They usually result 
in conservative estimations. Margolis and Vogt [29] used the recursive method of 
Zubov for 2nd-order systems. The approach is restricted to systems with nonlin- 
earities that possess Taylor series expansion and suffers from non-uniform conver- 
gence. Davison and Kurak [12] took the constrained minimization approach to fit 
a hyperellipse in the domain of attraction. Vannelli and Vidyasagar [51] used the 
concept of maximal Liapunov function to obtain a new partial differential equa- 
tion. The partial differential equation provides a basis for an iterative procedure 
to compute the domains of attraction. 

The complexity of these methods scales with the order of the systems and the 
order of the polynomials assumed for the Liapunov function. Manually computing 
the coefficients of the polynomials is very tedious. MACSYMA-like symbolic pack- 
ages employing computer algebra techniques can automatically determine these 
coefficients [42]. Much work remains in the area of applying symbolic methods to 
the construction of the Liapunov functions. 

6.3 Phase-Space Nonlinear Control 

The gain scheduling and cellular control approaches use, to some extent, the phase- 
space representation of control systems to design nonlinear control laws. Like Phase 
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Space Navigator, they discretize a phase space into regions; but their decompo- 
sitions are different and therefore bring to bear a different set of computational 
techniques. 

In gain scheduling [45], the system is linearized around a priori selected oper- 
ating points. A linear controller is designed for each operating point. Interpolation 
is used in regions between the operating points. Gain scheduling is an open-loop 
adaptation where parameters of a regulator change in a prespecified way. Autopi- 
lots and chemical process control are instances of gain scheduling. The method 
requires designers to decompose a phase space into locally linear regions. Phase 
Space Navigator differs from gain scheduling in that Phase Space Navigator auto- 
matically decomposes a phase space, even into regions that are global and nonlin- 
ear, and then synthesizes global reference trajectories using the knowledge of the 
phase-space decomposition. 

We have compared the method of the Phase Space Navigator with the dynamic 
programming approach in Section 4.3.4. The cellular control method of Hsu [20] 
shares very much the same disadvantages as dynamic programming. The method 
is based on the cell state space concept used in the cell-to-cell mapping. It is similar 
to dynamic programming in discretizing the state space, but differs in that it also 
partitions the cost function space. The control law is synthesized while certain 
cost functions are minimized. The control value at each cell state is programmed 
into a table. Run-time control is to perform a table-lookup. This method requires 
an extensive search for control paths at very fine grain. The number of cells that 
the cellular method needs to visit could be very large even in a moderately-high- 
dimensional state space. The Phase Space Navigator, however, decomposes the 
phase space into a manageable collection of nonuniform subregions, the flow pipes, 
and searches for global paths in this collection. 

Chaotic dynamical systems have extremely rich behaviors that can be tapped to 
design high-quality control systems. Recently, there have been several approaches 
to the controlling of chaos. Bradley developed a program called Perfect Moment 
that designs control paths for chaotic systems in phase space [8]. The program 
identifies phase-space chaotic features like strange attractors and searches for useful 
trajectory segments in a gridded phase space discretized into cells. Both Phase 
Space Navigator and Perfect Moment construct segmented control paths between 
the origin and destination states. Perfect Moment can be very useful in designing 
control paths for chaotic systems where the local structures of phase spaces are 
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extremely complicated. On the other hand, Perfect Moment could not characterize 
the gross geometric features of a chaotic attractor, say the envelope, and use the 
global description of the attractor to guide the search; it relies entirely on the 
mapping of the grid of cells. Phase Space Navigator assumes there is no chaos. 
It searches in a space of flow pipes, as opposed to a grid of cells. It explicitly 
models the global geometries of phase space and trajectory flows. Perfect Moment's 
gridded search could be combined with Phase Space Navigator's global phase-space 
modeling to enhance the search for control trajectories. 

The control design method developed in this thesis is not meant to replace other 
nonlinear control design techniques. To the contrary, the method complements 
the other techniques. The potential of the method in solving real problems can 
be fully realized when the method works together with other techniques. For 
example, feedback linearization is a well-developed nonlinear design method [45]. 
The method introduces a state feedback to a nonlinear system to cancel out the 
nonlinearity, thereby transforming the system into a linear one. The standard 
linear control design can be performed on the linear model. This linear design is 
then transformed back to work for the nonlinear system. Feedback linearization 
requires that the control input is continuous and the system has no zero dynamics. 
In cases when these conditions are met and a feedback linearizing control can be 
found, the feedback linearization could be an attractive alternative. On the other 
hand, Phase Space Navigator would be superior when, for example, the control 
takes discrete values as in switching-mode control. 

6.4 Intelligent Control 

The area of intelligent control, sometimes called knowledge-based control, has been 
very active recently. The knowledge-based approaches include direct expert control 
and indirect supervisory control [27]. The direct expert control encodes simple 
control laws into a rule base and invokes them at the run-time. The indirect 
approach influences the process by switching on and off controllers and updating 
control laws. For example, extreme alarm supervision switches on and off adaptors 
by monitoring the input signals and state variables. An intelligent control system 
usually has a hierarchical structure. The hierarchy ranks from high-level decision 
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making to low-level PID 1 controls. 

Phase Space Navigator can be regarded as one kind of intelligent control, if one 
wishes. It encodes deep domain knowledge of phase space and dynamics in the 
form of productions rules, algorithms, and other data structures. Different from 
the hierarchical structure in intelligent control, Phase Space Navigator consists of a 
global path planner and a local trajectory tracker; the two modules are integrated 
in a much tighter loop. This tight loop allows control decisions to be revised 
more naturally and less costly. For example, when the local tracker cannot correct 
trajectory deviation with local means, it calls upon the global planner to synthesize 
a new reference trajectory. An "intelligent controller" would have to pass the 
request and transform the decision back and forth through the control hierarchy. 



^ID stands for proportion, integral, and derivative. 
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Chapter 7 



Conclusions 



7.1 Thesis Revisited 

This thesis has developed a flow-pipe based phase-space method for designing non- 
linear controllers and a qualitative representation for encoding dynamics. It has 
formulated the task of control design and analysis as a computational one in which 
computations and reasoning about dynamics are essential. It has constructed a 
Workbench comprising the programs MAPS and Phase Space Navigator to as- 
sist control engineers and to automate a significant portion of the control design 
task. The programs have been applied to the design of a nonlinear controller for a 
magnetic levitation vehicle. 

MAPS analyzes the phase-space structure of a system in terms of geometric 
features. It produces a high-level description of the phase space that can be used 
to focus and prune the search for control paths. Because of its human accessible as- 
pect, the description provides meaningful information to control engineers in their 
design process. Phase Space Navigator automatically synthesizes control reference 
trajectories using the phase-space description as a "map". The "map" models the 
phase space with trajectory flow pipes and provides an efficient representation for 
the search. We described the algorithms for finite- valued control synthesis and 
mentioned caveats in generalizing them to continuous, multiple control parameter 
spaces in Section 4.3.4. When the computations are tractable, the method has 
great leverage for grossly nonlinear systems. 

The Workbench embodies domain knowledge from control engineering and dy- 
namical systems theory. It understands concepts like asymptotic and transient 
behaviors, stability regions, reachable sets, convergence, overshooting, etc. The 
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Workbench organizes the modules of analysis, design, and graphics presentation 
around the control task and encourages incremental changes to the Workbench, for 
example, incorporating programs tailored to particular applications and encoding 
knowledge of specific domains. Suppose we want to use the Workbench to design 
electric power control systems. In addition to their differential equation models, the 
power systems have special structural properties and model formulations that can 
be exploited to reduce design complexity and improve design quality. Special pro- 
gram fragments exploiting these properties can be integrated into the Workbench. 
The integration is made possible and easier by the underlying Scheme implemen- 
tation as a substrate. The Scheme programming language facilitates composition 
and abstraction of procedures. 

7.2 Problems and Future Work 

We examine the situations in which the Workbench does not work, analyze the 
reasons for its breakdown, and suggest further enhancements to the Workbench 1 
and future areas where much exciting work remains to be done. 

MAPS performs stability and flow analysis for a continuous dynamical system 
in a bounding box. It assumes that the system contains no chaotic attractors and 
fractal boundaries. MAPS can also analyze steady-state trajectories for discontin- 
uous systems; in fact, it has successfully analyzed the limit-cycle behaviors for a 
switching power regulator, a series resonant converter with clamped tank capacitor 
voltage, under various operating conditions [25]. 

To extend the stability and flow analysis of MAPS to work for systems with 
discontinuities is difficult. One reason is that the theoretical basis of the sta- 
bility analysis breaks down. Theoretical results for stability characterizations of 
discontinuous systems are scarce [16]. Computational simulation and modeling 
techniques, combined with ad hoc treatments for limit cycles, could provide an 
attractive alternative for analyzing discontinuous systems. 

Chaotic attractors are difficult to characterize. The boundaries and the attrac- 
tors themselves can be fractal. Their fractal boundaries are typically caused by the 
heteroclinic intersections of stable and unstable manifolds. The fine, convoluted 
structures pose great problems for numerical algorithms attempting to determine 



1 Chapter 3 has discussed several extensions to the MAPS program. 
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them. The stability and robustness of the computation is the main issue in this 
case. 

The Workbench could be extended to characterize gross geometric features of 
chaotic attractors. It could augment purely numerical measures with some invari- 
ant, symbolic characterizations. Many chaotic attractors, for example the Rossler 
band [5], are constructed with a stretching- and-folding operation. The Workbench 
could compute the Birkhoff signature encoding the heteroclinic intersections of sta- 
ble and unstable manifolds. This symbolic sequence plus the stretching-and-folding 
operation can help understand the fractal structures. The Workbench could also 
calibrate dimensions of the attractors and determine Liapunov exponents. The 
combination of geometric descriptions, symbolic constraints, and numerical mea- 
sures can better characterize the fine structures of chaotic attractors. 

The Workbench has adopted the Euclidean metric in n-dimensional Euclidean 
spaces. However, a different metric, the one on surfaces of manifolds rather than 
on their embedding spaces, could also be used. This new metric would be able 
to more accurately model the folding of a manifold. The information about the 
geometric properties of the manifold can be used to improve the geometric mod- 
eling algorithm discussed in Chapter 3. In addition to the metric information, the 
Workbench could use certain topological information about the space in planning 
control. We have implemented a program for computing the homology group of a 
triangulated manifold based on a reduction algorithm [33]. It will be interesting to 
see how we could use this topological property of phase portraits to design control 
paths. The Workbench could also be extended to work for other types of spaces 
such as projective spaces and tori. 

The design capability of Phase Space Navigator can be enhanced with other 
optimization techniques. For example, the variational method can be used in 
optimizing the local control trajectories after the global path is established. The 
cost of the optimization at the stage of local trajectory generation is expected to 
be low. 

Future research should address the robustness of the phase-space control de- 
sign regarding for example model order uncertainties or unmodeled dynamics. The 
technique of "thickened trajectory" discussed in Section 4.3 models only measure- 
ment uncertainties or noise. The structural uncertainties of the model could change 
the topological structure of a phase space. The effect of such uncertainties on the 
phase-space reference trajectories remains to be quantified. The thesis has not 
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addressed model building for dynamical systems. The Workbench could interface 
with a modeling program that allows the Workbench to deal with a wider class of 
systems. The modeling program should formulate the most appropriate model for 
the Workbench. 

The current implementation of the Workbench is expected to work well for mod- 
erately low-dimensional systems. However, "the curse of dimensionality" makes 
the computations intractable for high-dimensional systems. This is somewhat a 
paradox: high-dimensional, highly-nonlinear systems exhibit richer flow patterns 
in phase space and provide more design opportunities for the phase-space method. 
However, such advantage is at the expense of more computations. To extend the 
applicable range of the method, the efficiency of the geometric algorithms needs 
to be greatly improved. 

7.3 Broad Implications 

The work reported in this thesis has broad implications: 

The foremost is its demonstration that controllers for complex systems can be 
automatically designed. The work showed that programs can visualize qualitative 
behaviors in phase space; programs can plan control trajectories and steer systems; 
and the phase-space geometric modeling provides a "map" for navigating systems 
in phase space. 

Second, novel computational representation and reasoning mechanisms can be 
developed in the context of automating challenging engineering tasks. The dynam- 
ics of nonlinear systems is difficult to describe and manipulate. The qualitative 
representation developed in this thesis provides a way to computationally describe 
the qualitative aspects of the dynamics. With this representation, the difficult con- 
trol design is translated into a computational task: the flow-pipe based mechanism 
manipulates a system's natural dynamics and synthesizes the desired dynamics for 
the system. 

Third, the domain knowledge and techniques from symbolic, numerical, and 
geometric computing have been proven essential. In order to fully exploit the dy- 
namics and to build programs to imitate human control designers, the Workbench 
uses whatever knowledge and techniques that are necessary: geometric theory of 
dynamical systems, control theory, and techniques from artificial intelligence, vi- 
sion, computational geometry, numerical analysis, and graph search. 
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Fourth, the complexities of the control design task necessitate the need for 
automatic modeling and analysis tools. Autonomous programs like those in the 
Workbench have extended the capabilities of the "eyes" and "hands" of control 
engineers in seeing and manipulating objects. They have enlarged the design space 
engineers can explore. In certain cases, programs can even outperform human 
designers. 

Fifth, the phase-space trajectory planning has the potential of becoming a pow- 
erful new paradigm for nonlinear control synthesis and a complementary alternative 
to other control design techniques. 

Finally, the potential applications of the Workbench are numerous: (1) Large- 
scale power systems are difficult to analyze and control. A power system needs 
to be analyzed in terms of its stability margin. It needs to be brought back to 
the original equilibrium after a power outage. The Workbench could be enhanced 
to provide geometric means for the stability analysis and possibly for the control 
trajectory design. (2) Vehicles like the maglev trains and automobiles usually use 
passive suspension to dampen disturbances caused by road irregularities and wind 
loads. It is conceivable that such disturbances could be compensated with active 
actuation. This technology is called active suspension. In the maglev example, it is 
possible to characterize the guideway irregularities; thus active suspension could be 
an attractive candidate. The Workbench can help in designing such a secondary 
suspension for the maglev vehicle that will complement the primary suspension 
(magnetic levitation) designed in Chapter 5. (3) Large flexible space structures 
use thin beams as building blocks in structures like the truss. These beams work 
in nonlinear regions. The Workbench can be used to study the buckling of the 
beams and explore control schemes to strengthen them. There are many more 
areas where the Workbench will find itself useful. 

7.4 Towards A Language for Computational Control De- 
sign 

We outline a computational language with which one can describe aspects of the 
control design. Although the language is in a very preliminary form, elucidating 
basic constructs and properties of the language helps capture the informal design 
procedures and serves as a stepping stone for formalizing the knowledge of control 
engineering. 
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Professional control engineers make use of implicit working knowledge such 
as the phase-plane method and work at different design levels towards the goal. 
By expressing the informal design knowledge explicitly in computational terms, 
we are able to build programs to manipulate the knowledge, mechanize certain 
laborious tasks, communicate the knowledge more effectively to other designers, 
and improve education in control engineering. In addition, the language provides 
means for building design abstractions, for capturing design rationales, and for 
addressing computational complexities of the design. 

The design language consists of a set of primitives and operations on the prim- 
itives. The primitives are phase portraits segmented into flow pipes. We call these 
segmented portraits the flow-pipe portraits. The operations combine these flow- 
pipe portraits to form new ones. An operator also takes conditional choice as an 
argument. In composing two flow-pipe portraits, the operator intersects the flow 
pipes in the portraits and glues them together at intersections. The choice infor- 
mation specifies for flow pipes which way to go at each intersection. The semantics 
of the language is dictated by the underlying dynamics of the system and the con- 
trol task. For example, an intersection of two flow pipes is specified as a geometric 
intersection of their polyhedral approximations. 

We express the phase-space control design in a high-level form with this lan- 
guage: the control design is the composition of flow-pipe portraits. The opti- 
mization of the design is equivalent to the minimization of certain measure in the 
composition process. When an optimization criterion is specified, the composition 
process constructs the flow-pipe portrait containing the optimal control path. 

7.5 The Big Picture 

The domain of control synthesis is interesting to study for several reasons. First, 
the computational codification of engineering knowledge such as the knowledge 
of control design enables us to distill and communicate principles of engineering 
design in a more explicit way. Second, practical control systems are extremely 
useful but are often difficult to design. Intelligent design aids for control engineers 
like the Control Engineer's Workbench can improve the design quality and increase 
the productivity. 

I am convinced that high-performance computational techniques have great 
potentials in revolutionizing traditional engineering computation and design, and 
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research investigating challenging real-world problems serves as a driving force for 
the development of computer science and artificial intelligence. This thesis serves 
as an instance for demonstrating my conviction. 

This work has primarily focused on the automatic synthesis of global nonlin- 
ear controllers and has not addressed the important issues of modeling, sensing, 
estimation, and high-quality linear controller design. The exploration and de- 
velopment of a collection of new methods tackling these issues in all dimensions 
can revolutionize the synthesis and analysis of high-performance nonlinear control 
systems. 

The research described in this thesis actively exploits new synergies among 
artificial intelligence, computer science, applied mathematics, dynamical system 
theory, and control engineering and develops appropriate computational techniques 
to tackle real problems in engineering. The Control Engineer's Workbench is a 
prototype of a new class of intelligent computational tools that combine numerical 
and symbolic computations with AI reasoning techniques and automatically model, 
analyze, and design complex physical systems. These tools, when used as human 
aids, will greatly amplify and extend our capabilities in seeing and manipulating 
the world around us. 
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Appendix A 

Maglev Application: Equilibria 



The equilibria for V { = 140 and V { = 300: 



; Vi = 140 
((#(140. 0. 200.) 
(-17. 004+22. 963i 
-17. 004-22. 963i 
24.007) 
(#(.23604 .97174 0) 
#(.51331 -.55588 .65384) 
#(.30157 .73255 .61027)) 
saddle)) 



equilibrium point #(x y z) 
eigenvalue 1 
eigenvalue 2 
eigenvalue 3 
eigenvector 1 
eigenvector 2 
eigenvector 3 
stability type 



; Vi = 300 
((#(300. 0. 428.57) 
(-21. 411+21. 394i 
-21. 411-21. 394i 
21.394) 
(#(.54819 .83636 0) 
#(.43085 -.65947 .616) 
#(.23229 .71052 .66423)) 
saddle)) 



equilibrium point #(x y z) 
eigenvalue 1 
eigenvalue 2 
eigenvalue 3 
eigenvector 1 
eigenvector 2 
eigenvector 3 
stability type 
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Appendix B 

Maglev Application: Synthesized 
Reference Trajectories 



This appendix contains the initial, switching, and final points of the synthesized 
reference trajectories starting from four different initial states, as shown in Fig- 
ure 5.6. Each point in the print-out is a vector of (t,x,y,z). 

; trajectory 1 

; init-pt: #(0. 140. 0. 220.) 

((time 0.) (switching-state #(140. 0. 220.)) 

(control 300.)) 
((time .013441) (switching-state #(162.98 4.4416 221.42)) 

(control 140.)) 
((time .12704) (switching-state #(139.39 .078872 201.75)) 

(control *local-control*)) 

; trajectory 2 

; init-pt: #(0. 140. 0. 280.) 

((time 0.) (switching-state #(140. 0. 280.)) 

(control 300.)) 
((time .10156) (switching-state #(276. -74.214 338.79)) 

(control 140.)) 
((time .32041) (switching-state #(139.26 -.19212 196.51)) 

(control *local-control*)) 
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; trajectory 3 

; init-pt: #(0. 140. 0. 290.) 

((time 0.) (switching-state #(140. 0. 290.)) 

(control 300.)) 
((time .19536) (switching-state #(264.69 -105.03 340.68)) 

(control 140.)) 
((time .42865) (switching-state #(142.4 3.2661 202.69)) 

(control *local-control*)) 

; trajectory 4 

; init-pt: #(0. 140. 0. 300.) 

#(0. 140. 0. 300.) is outside controllable region 
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